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Abstract 

This  research  presents  an  algorithm  that  improves  the  ability  to  view  objects 
using  an  electro-optical  imaging  system  with  at  least  one  polarization  sensitive  channel 
in  addition  to  the  primary  channel. 

Following  the  review  of  historical  methodologies  applicable  to  this  research  area, 
the  statistical  Cramer-Rao  lower  bound  (CRLB)  is  developed  for  a  two-channel  po- 
larimeter.  The  CRLB  is  developed  using  the  system’s  ability  to  resolve  two  point 
sources  in  the  presence  of  atmospheric  turbulence.  The  bounds  show  that  such  a 
polarimeter  has  an  advantage  over  previous  imaging  methods  at  smaller  separations. 

A  small  optical  laboratory  is  set  up  to  generate  a  set  of  calibrated  images  for  ver¬ 
ification  of  the  simulation  results  and  validation  of  algorithm  development.  Defocus 
is  the  aberration  chosen  for  algorithm  development  and  testing  clue  to  its  significant 
presence  when  imaging  through  turbulence  and  its  ease  of  production  in  the  labora¬ 
tory.  An  innovative  algorithm  for  detection  and  estimation  of  the  defocus  aberration 
present  in  an  image  is  also  developed. 

Using  a  known  defocus  aberration,  an  iterative  polarimeter  deconvolution  algo¬ 
rithm  is  developed  using  a  generalized  expectation-maximization  (GEM)  model  that 
produces  results  as  predicted  by  the  CRLB  results.  Using  an  example  bar  target 
set  with  a  degree  of  polarization  of  one,  the  polarimeter  deconvolution  algorithm  can 
resolve  the  two  bars  down  to  half  the  bar  separation  as  the  Richardson-Lucy  (RL) 
algorithm  can  do.  In  addition,  a  fidelity  metric  is  used  that  shows  the  polarimeter 
deconvolution  algorithm  deconvolves  simulated  targets  with  approximately  half  of  the 
error  present  in  objects  deconvolved  using  the  RL  algorithm. 

The  polarimeter  deconvolution  algorithm  is  extended  to  an  iterative  polarime¬ 
ter  multiframe  blind  deconvolution  (PMFBD)  algorithm  with  an  unknown  aberration. 
Using  both  simulated  and  laboratory  images,  the  results  of  the  new  PMFBD  algo- 


rithm  clearly  outperforms  an  RL-based  MFBD  algorithm.  The  convergence  rate  is 
significantly  faster  with  better  fidelity  of  reproduction  of  the  targets. 

This  research  successfully  developed  an  algorithm  that  uses  polarization  data  in 
conjunction  with  standard  imaging  to  improve  the  spatial  resolution  of  deconvolved 
objects  with  faster  convergence  rates.  Clearly,  leveraging  polarization  data  in  electro- 
optical  imaging  systems  has  the  potential  to  significantly  improve  the  ability  to  resolve 
objects  and,  thus,  improve  Space  Situation  Awareness. 
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P OLARIMETER  BLIND  DECONVOLUTION 

Using  Image  Diversity 

I.  Introduction 

As  telescopes  have  grown  larger  and  larger  the  effect  of  atmospheric  turbulence 
on  the  astronomical  images  has  increased.  If  the  distortion  introduced  by  the 
atmosphere  on  the  incoming  light  can  be  determined,  it  can  be  removed  from  the  im¬ 
age.  This  research  addresses  an  approach  for  improving  image  resolution  through  the 
use  of  a  polarimeter  and  the  diversity  in  the  images  it  produces.  A  polarimeter  pro¬ 
duces  multiple  simultaneous  observations  of  the  target  that  are  related.  The  difference 
between  the  images  is  due  to  the  polarized  light  present  in  the  different  channels.  The 
polarized  light  is  created  through  the  interaction  of  light  and  materials.  Additionally 
these  observations  are  distorted  by  the  same  atmospheric  turbulence. 

The  traditional  blind  deconvolution  problem  is  ill-posed  in  that  both  the  system 
response  and  the  input  image  are  to  be  determined  from  a  single  composite  observa¬ 
tion.  Historical  deconvolution  methods  use  iterative  approaches  to  arrive  at  a  solution 
to  the  problem.  The  ability  to  improve  the  solution  is  the  main  motivation  for  this 
research 

1.1  Space  Situational  Awareness 

Any  improvement  in  the  ability  to  resolve  objects  viewed  through  atmospheric 
turbulence  improves  Space  Situational  Awareness  (SSA).  Its  importance  is  represented 
in  the  following  two  quotes.  The  first  is  an  excerpt  from  a  speech  by  Lt  Gen  Klotz, 
Vice  Commander,  AFSC  in  2006. 

In  the  face  of  potentially  emerging  threats,  our  fundamental  U.S.  position 
is  straightforward  and  our  new  National  Space  Policy  lays  it  out  very 
well,  and  I  commend  it  to  your  reading.  Let  me  try  to  summarize  it  in 
this  way:  We  must  ensure  access  to  space,  as  it  is  critical  to  our  Nation’s 
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full  range  of  capabilities  ...  capabilities  that  transform  the  battlefield  and 
empower  the  global  economy.  Indeed,  the  United  States  is  committed  to 
the  free  use  and  access  to  space  by  all  nations  for  peaceful  purposes.  The 
reality  is,  however,  we  cannot  assume  all  nations  will  pursue  only  peaceful 
purposes  ...  just  the  same  way  we  can  not  assume  all  nations  will  respect 
the  domains  of  land,  air  or  sea.  In  light  of  this  reality,  we  must  protect 
our  key  capabilities.  They  are  a  vital  National  interest  and  one  of  our  top 
priorities. 

The  first  critical  component  of  that  protection  is  space  situational  aware¬ 
ness.  Currently,  we  in  Air  Force  Space  Command  catalog  approximately 
14,000  objects  in  space  ...  but  it’s  just  that:  a  catalog.  While  we  can 
detect  objects  as  small  as  a  baseball  in  low  Earth  orbit  ...  or  as  small  as 
a  basketball  in  geosynchronous  orbit,  we  need  more  than  just  the  location 
and  general  direction  of  an  object.  Our  aim  is  to  move  beyond  cataloging 
to  understanding  what  is  “up  there.”  [23] 

The  second  is  an  excerpt  from  the  June  2006  Air  Force  Magazine  [44]: 

...The  new  Quadrennial  Defense  Review  (QDR),  released  in  February,  sim¬ 
ply  noted  that  Washington  must  have  “unfettered,  reliable,  and  secure” 
access  to  its  space  assets,  assured,  for  now,  by  “improving  space  situa¬ 
tional  awareness  and  protection,  and  through  other  space  control  mea¬ 
sures.”  The  Air  Force  is  taking  its  cue  from  the  QDR,  focusing  most  of 
its  nonclassified  efforts  at  space  superiority  on  systems  that  will  broadly 
enhance  its  knowledge  of  what’s  in  orbit,  as  well  as  its  ability  to  know  if 
American  space  systems  are  under  attack.  ...  “We  have  to  know  what’s 
up  there,”  said  Gen.  T.  Michael  Moseley,  Air  Force  Chief  of  Staff.  “We 
have  to  continually  modernize  the  early  warning  systems  to  know  what  is 
up  there,  what  has  been  added,  what  are  the  orbital  paths,  and  what  are 
the  opportunities  to  see.”  This  is  what  the  United  States  must  do  to  avoid 
“a  Pearl  Harbor  in  space,”  Moseley  observed.  The  emphasis  remains  on 
space  situational  awareness... 


The  main  motivation  for  this  research  is  in  improving  that  SSA  capability. 


1.2  Problem  and  Definitions 

1.2.1  Geometry  of  Scenario.  The  scenario  of  interest  involves  the  incoherent 
illumination  of  satellites  by  sunlight  which  is  then  reflected  by  the  satellite  surfaces 
and  viewed  by  a  telescope.  The  light  reflected  can  be  partially  polarized.  The  amount 
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of  polarization  is  determined  by  the  incident  angle  of  the  light  and  the  material  that 
reflects  the  light  as  described  in  more  detail  in  Section  2.1. 

1.2.2  Imaging  .  Imaging  systems,  such  as  telescopes,  are  used  to  form 
images  of  objects.  The  incoming  light  from  an  object  is  transformed  by  a  lens  to  form 
an  image  of  the  object. 


Object  — >  Imaging  System  —>■  Image 


The  intensity  distribution  of  the  object  is  convolved  with  the  point  spread  func¬ 
tion  (psf)  to  form  the  image.  Using  Fourier  transforms  to  move  from  the  spatial 
domain  into  the  spectral  domain  changes  the  convolution  into  a  multiplication.  The 


Fourier  transform  of  the  psf  is  the  optical  transfer  function  (OTF).  These  relationships 


can  be  seen  in  Figure  1.1. 
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Figure  1.1:  Fields  and  Relationships. 
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1.2. 2.1  Diffraction-limited  Image.  In  a  vacuum,  the  image  formed  is 
limited  only  by  diffraction  of  the  light  as  it  passes  through  an  ideal  imaging  system. 
Figure  1.2(a)  shows  the  OTF  of  a  circular  aperture  and  Figure  1.2(b)  shows  the 
resulting  image  of  a  point  source. 


OTF  OTF  Cross  Section 


(a)  OTF  of  Circular  Aperture. 


Detector  Image 
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Detector  Image  Cross  Section 


(b)  Point  Source  Image. 


Figure  1.2:  OTF  and  Image  with  Cross  Sections. 


1.2. 2. 2  Atmospheric  Turbulence.  Atmospheric  turbulence  spatially 
distorts  the  wavefront  as  light  passes  through  it  and  causes  blurring  of  images  in  an 
imaging  system.  Figure  1.3  shows  an  example  of  distortion  caused  by  a  randomly 
generated  phase  screen  using  10  Zernike  coefficients.  The  upper  left  image  is  the 
resulting  phase  screen.  The  upper  right  image  is  the  OTF  created  using  the  phase 
screen.  The  lower  set  of  images  are  of  a  point  source  using  the  OTF  and  its  cross 
section. 


1.2. 2. 3  Adaptive  Optics.  Adaptive  optics  (AO)  systems  are  used  to 
reduce  the  effects  of  atmospheric  turbulence.  Using  deformable  mirrors,  an  AO  system 
can  detect  and  significantly  reduce  lower  order  terms  of  the  spatial  distortion.  Two 
of  the  lowest  order  terms,  tip  and  tilt,  account  for  approximately  80  percent  of  the 
magnitude  of  the  distortion  and  can  be  significantly  reduced  by  a  movable  mirror. 
The  next  several  terms  are  typically  reduced  by  a  separate  deformable  mirror.  The 
remaining  terms  require  post-processing  to  reduce  their  impact  on  the  image.  Figure 
1.4  shows  an  example  of  the  problem.  The  left  image  is  a  computer  generated  image  of 
the  Hubble  Space  Telescope  (HST).  The  right  image  is  from  a  3.7m  Advanced  Electro- 
Optical  System  (AEOS)  telescope  at  the  Air  Force  Maui  Optical  and  Supercomputing 
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(a)  Phase  Screen 


(b)  OTF  with  phase  distortion 


(c)  Det.  Image  with  Distortion  (d)  Det.  Image  with  Dist.  Cross-section 


Figure  1.3:  OTF  and  Detector  Image  with  10  Zernikes. 

Site  (AMOS)  with  adaptive  optics  on  and  tracking  properly.  As  can  clearly  be  seen, 
or  not  as  the  case  may  be,  there  is  a  significant  amount  of  atmospheric  distortion 
remaining  in  the  image. 

1.2.3  Convolution.  The  relationships  between  the  observed  image  and  the 
actual  object  can  be  defined  mathematically.  In  imaging,  the  observed  image  i(x,  y )  is 
the  two-dimensional  convolution  of  the  true  image  of  the  object  o(z ,  w)  with  a  linear 
shift-invariant  blur,  the  psf,  h(x,y ): 


i(x,y)  =  h(x,  y)  ©  o(x,  y), 


(1.1) 

(1.2) 
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(a)  Computer  Model  of  the  (b)  Example  Diffraction  Lim-  (c)  Real  Image  of  the  Hubble 
Hubble  Space  Telescope.  ited  Image.  Space  Telescope  from  3.7m  AEOS 

Telescope. 


Figure  1.4:  Hubble  Space  Telescope  Images. 


where  ©  is  the  two-dimensional  convolution  operator. 


The  transfer  function  between  the  Fourier  transform  of  the  object  viewed  and 
the  Fourier  transform  of  the  image  captured  by  an  imaging  system  is  the  OTF  [15]. 
The  relationships  between  the  object,  its  image,  and  the  Fourier  domain  are  shown 
in  Figure  1.5.  Using  these  relationships,  it  can  be  seen  that  the  image  is  the  inverse 
Fourier  transform  of  the  Fourier  transform  of  the  object  multiplied  by  the  OTF. 


hbject  ©  PSF 
Sobject  X  OTF 


limaqe 

T 

Simage 


Figure  1.5:  Image,  Object,  Fourier  Domain  Relationships. 


lobject  and  I  image  are  the  intensity  distributions  of  the  object  and  image,  respec¬ 
tively.  S0bject  and  Simage  are  their  corresponding  Fourier  transforms  in  the  spectral 
domain. 


1.2.4  Deconvolution.  Deconvolution  is  removing  the  distortion  from  the 
image  by  dividing  the  Fourier  transform  of  the  image  by  the  Fourier  transform  of  the 
distortion  which  is  the  optical  transfer  function.  This  ideal  process  is  shown  in  Figure 
1.6.  If  the  OTF  is  known  and  there  is  no  noise  then  the  deconvolution  is  nearly 
perfect.  Only  the  error  from  the  Fourier  transforms  attributable  to  finite  machine 


6 


precision  remains.  The  problem  is  that  in  many  cases  the  OTF  that  created  an  image 
can  never  be  known  but  only  estimated  from  the  intensity  data.  The  reason  that 
iterative  approaches  are  desireablc  is  the  presence  of  zeros  in  the  OTF  make  division 
by  a  real  OTF  impossible. 


L  image 


IT 


I object 

T 


S, 


image 


/OTF  =  S, 


object 


Figure  1.6:  Deconvolution  Process. 


Several  blind  deconvolution  methods  have  been  developed  to  estimate  the  OTF 
from  images  that  were  created  by  unknown  psfs  and  are  discussed  in  Section  2.2. 


1.2.5  Polarimeter  Imaging.  A  very  basic  polarimeter  is  constructed  using  a 
polarizing  beam  splitter  (PBS).  If  incident  light  with  a  vertical  polarization  is  trans¬ 
mitted  by  the  PBS  then  incident  light  with  a  horizontal  component  is  reflected.  This 
results  in  two  simultaneous  images  differing  only  by  the  polarization  components.  The 
light  from  the  object,  o(z,w),  is  a  composite  of  the  orthogonal  components,  Oh(z,w ) 
and  ov(z,w).  Therefore  o(z,w )  =  Oh(z,w)  +  ov(z,w).  Using  this  equation  and  ap¬ 
plying  Equation  1.1  to  a  polarimeter  results  in  the  corresponding  intensity  images, 
ih(z,w )  and  iv(z,w),  such  that 


and 


ih(x,  y)  =  h(x  -  u  y  ~  w)°h(zi  w)i 

z  w 

(1.3) 

iv(x,  y)  =  ^2^2h(x  -  z,y  -  w)ov(z ,  w). 

(1.4) 

Z  W 


The  diversity  between  the  two  images  is  created  by  the  polarizing  effects  of  the 
object  observed.  The  wavefront  for  both  images  is  assumed  to  be  identically  distorted 
by  the  atmospheric  turbulence  and,  thus,  the  same  psf  h(x,  y )  is  used.  This  reduces 
the  number  of  unknowns  between  the  two  equations  and  is,  thus,  an  advantage  over 
traditional  blind  deconvolution.  Developing  a  useful  relationship  between  Oh(z,w ) 
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and  ov(z,w)  would  further  reduce  the  number  of  unknowns  to  be  estimated  from  the 
measured  data. 

1.3  Layout  of  the  Dissertation  document 

The  goal  of  this  research  is  to  understand  the  issues  involved  with  viewing 
objects  with  a  polarimeter  and  develop  a  method  to  improve  spatial  resolution  with 
the  additional  polarity  data  that  it  provides.  Three  objectives  are  defined  to  achieve 
this  goal  followed  by  the  results  of  the  algorithm  development. 

1.3.1  Polarimeter  Spatial  Resolution  Bounds.  The  first  objective  is  an 
investigation  of  the  Cramer-Rao  Lower  Bound  (CRLB)  for  spatial  resolution  of  a  two- 
channel  polarimeter.  The  CRLB  is  used  to  determine  if  the  additional  polarization 
information  provides  any  grounds  for  improvement  in  spatial  resolution  over  standard 
imaging  techniques.  The  bounds  are  built  in  several  steps.  Starting  with  development 
of  the  CRLB  of  a  standard  imaging  system  without  atmospheric  effects  [43]  and 
building  to  the  CRLB  of  a  two-channel  polarimeter  incorporating  atmospheric  effects 
[42] .  The  CRLB  development  is  presented  in  Chapter  111  along  with  a  comparison  of 
the  polarimeter  CRLB  to  standard  imaging  systems.  Understanding  the  conditions 
in  which  the  polarimeter  can  outperform  standard  imaging  systems  is  the  desired 
outcome  for  the  objective. 

1.3.2  Laboratory  and  Simulation  Data.  The  second  objective  is  to  pro¬ 
duce  calibrated  data  for  testing  any  developed  algorithms.  The  laboratory  setup  used 
allows  the  controlled  introduction  of  focus  aberrations  in  the  optical  path.  An  innova¬ 
tive  maximum -a-priori  (MAP)  estimation  approach  is  developed  in  Section  4.2  that 
estimates  the  focus  aberration  present  from  a  single  image.  The  resulting  images  are 
then  used  to  set  the  statistical  properties  of  the  Matlab®  simulation  and  validate  the 
performance  of  the  deconvolution  algorithms.  The  images  used  in  development  and 
testing  of  the  polarimeter  deconvolution  algorithms  are  discussed  in  Chapter  IV. 


1.3.3  Polarimeter  Deconvolution  Algorithms.  The  third  objective  is  to  de¬ 
velop  a  method  that  uses  the  polarization  data  to  improve  spatial  resolution.  Initial 
development  is  under  the  assumption  of  a  known  psf  with  the  results  compared  to 
a  traditional  deconvolution  algorithm  using  a  metric  based  on  the  fidelity  of  the  de¬ 
convolved  object.  The  algorithm  is  then  extended  to  a  polarimeter  multiframe  blind 
deconvolution  approach  similar  to  Schulz  [36]  but  incorporating  polarization  informa¬ 
tion.  The  polarimeter  deconvolution  algorithm  development  is  discussed  in  detail  in 
Chapter  V. 

1.3.4  Results.  Using  both  Matlab®  simulation  and  calibrated  laboratory 
images,  the  results  of  the  algorithm  development  are  detailed  in  Chapter  VI.  The 
significant  improvements  in  fidelity  of  object  deconvolution  and  convergence  rates 
when  compared  to  traditional  methodologies  are  presented.  Results  for  polarimeter 
deconvolution  are  presented  for  both  known  and  unknown  aberrations. 

1.3.5  Conclusions  and  Recommendations.  The  conclusions  drawn  from  the 
results  of  the  algorithm  development  are  presented  in  Chapter  VII  along  with  recom¬ 
mendations  for  follow-on  research. 
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II.  Historical  Background 


This  chapter  describes  key  concepts  involved  in  developing  a  method  to  improve 
imaging  using  polarimeters.  Starting  with  a  brief  description  of  a  polarimeter 
and  the  properties  of  materials  that  produce  polarized  light.  Section  2.2  follows 
the  development  of  historical  methods  that  led  to  blind  deconvolution  of  imaging 
through  turbulence.  The  final  section  briefly  describes  some  acceleration  techniques 
and  bounds  on  imaging. 

2.1  Polarimeters 

A  basic  polarimeter  is  formed  by  splitting  the  incoming  light  with  a  PBS  shown 
in  Figure  2.1.  The  PBS  splits  the  light  into  two  orthogonally  polarized  beams.  In 
this  example,  the  reflected  beam  is  arbitrarily  designated  to  be  vertically  polarized, 
Vp0i,  while  the  transmitted  beam  is  designated  to  be  horizontally  polarized,  Hpoi. 

Incoming  Light 


Vpol 


Hpoi 

Figure  2.1:  Basic  Polarimeter. 

2.1.1  Material  Reflectance  Models.  An  ideal  PBS  splits  unpolarized  light 
equally.  However,  when  unpolarized  light  is  reflected  from  certain  materials  it  becomes 
partially  polarized.  A  considerable  amount  of  data  is  available  on  the  polarization 
properties  of  materials  clue  to  reflectance.  Active  imaging  polarization  studies  char¬ 
acterize  properties  of  materials  by  varying  the  polarization  states  of  the  illumination 
and  analyzing  the  reflected  image.  Some  of  the  basic  principles  of  reflectance  from  a 
metal  are  discussed  in  the  following  section. 
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From  Hecht  [16]  the  amplitude  coefficient  for  reflection  perpendicular  to  the 
surface  is 


rii  cos  0,;  —  nt  cos  Qt 

r_ l  = - - - — 

n%  cos  0,j  +  nt  cos  Qt 

and  the  amplitude  coefficient  for  reflection  parallel  to  the  surface  is 


(2.1) 


nt  cos  (-)t  —  rii  cos  Qt 
nt  cos  (-)/  +  nt  cos  0  j 


(2.2) 


where  n  is  the  index  of  refraction  and  0,  is  the  angle  of  incidence  and  Qt  is  the  angle 
of  transmission.  The  subscripts  i  and  t  indicate  the  incident  or  transmissive  media  at 
the  reflective  interface. 

At  normal  incidence  0*  =  0,  and  ry  and  r±  simplify  to 


rk  ~  nt  ,  nt  -  n% 

r l  = -  and  r  n  = - . 

rii  +  nt  rii  +  nt 


(2.3) 


For  reflectance  R±  =  r]_  and  R\\  =  rjj.  Squaring  Equation  2.3  results  in  =  R\\  for 
normally  incident  light.  This  is  due  to  the  homogeneous  surface  of  the  metal. 


2. 1.1.1  Plasma  Frequency.  From  Hecht  [16]  “Free  electrons  and  pos¬ 
itive  ions  within  a  metal  may  be  thought  of  as  a  plasma  whose  density  oscillates  at 
a  natural  frequency  lop,  the  plasma  frequency”  The  plasma  frequency  is  related  to 
the  plasma  wavelength  through  the  equation  ujp  =  2i:c/\p  where  \p  is  the  plasma 
frequency  and  c  is  the  speed  of  light  in  a  vacuum.  The  dispersion  equation  using  u>p 
becomes 

nV)  =  1  -  K/^)2.  (2.4) 

The  following  can  be  observed  from  Equation  2.4: 

•  When  u  =  ujp  then  n(co)  =  0. 

•  When  tv  >  (jjp  then  n{eo)  is  real  and  the  metal  is  transparent.  (See  left  image  in 
Figure  2.2.) 
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When  ui  <  ujp  then  n(uj)  is  complex.  (See  right  image  in  Figure  2.2.) 


Figure  2.2:  Real  and  Imaginary  Components  of  n. 


Mandel  and  Wolf  [26]  state  that  the  relationship  between  the  refractive  index 
of  the  medium,  n(u),  and  the  dielectric  susceptibility,  n(u)  in  “a  medium  whose 
macroscopic  properties  do  not  change  in  time”  is 

n2(uj)  =  1  +  47 xn(u)  (2.5) 


where  a ;  is  the  frequency  of  interest.  (Note  the  different  notation  of  where  the  n  is 
located.) 

The  equation  from  Hecht  [16]  for  reflectance  is 


(nR-l)2  +  n2j 
(nR  +  l)2  +  n2j' 


(2.6) 
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Using  Equations  2.4  and  2.5  yields  the  following: 


n2{  u) 

=  1  +  47rn(u;) 

=  1  +  47t(1  -  (ujp/uj)2)^ 

(2.7) 

n{u) 

=  (1  +  4tt(1-(c np/u;)2)¥ 

(2.8) 

riR 

=  9Re(n(u;)) 

(2.9) 

ni 

=  5sm(n(u))) 

(2.10) 

Plugging  these  results  into  Equation  2.6  yields  the  reflectance  shown  in  Figure  2.3. 


Figure  2.3:  Reflectance  versus  Wavelength. 

If  one  were  to  adjust  the  plasma  frequency  to  approximately  0.31  /xm  and  add 
an  offset  of  n  =  0.35  (the  refractive  index  of  silver  [41])  to  the  refractive  index  given 
by  Mandel  and  Wolf  then  the  resultant  graph  (see  Figure  2.4)  resembles  the  plot  of 
the  reflectance  of  silver  as  reproduced  from  Figure  4.59  in  Hecht  [16]. 

These  models  relate  the  reflectance  of  different  materials  with  the  polarization 
states  produced  by  reflection  from  that  material.  This  supports  the  notion  that 
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Wavelength  (m) 

(c)  Silver 


Figure  2.4:  Predicted  Ur,  nj,  and  Reflectance  of  Silver. 


polarization  data  might  be  usable  in  identifying  materials  and/or  the  orientation  of 
the  surface  from  which  the  light  is  reflected. 


2.2  Blind  Deconvolution 

This  section  discusses  the  development  of  approaches  to  the  ill-posed  problem 
of  blind  deconvolution.  First  general  deconvolution  techniques  assuming  a  known 
psf  are  presented.  Then  approaches  for  blind  deconvolution  with  standard  imaging 
systems  are  presented,  along  with  a  few  means  of  accelerating  their  typically  lengthy 
execution  times.  Finally,  a  couple  of  papers  relatating  to  lower  bounds  of  resolution 
on  imaging  systems  are  presented. 

2.2.1  Single  Channel.  In  1972,  Gerchberg  and  Saxton  [13]  presented  an 
algorithm  “for  the  rapid  solution  of  the  phase  of  the  complete  wave  function  whose 
intensity  in  the  diffraction  and  imaging  planes  of  an  imaging  system  are  known.”  The 
algorithm  iterates  back  and  forth  between  the  two  planes  via  Fourier  transforms.  The 
problem  is  constrained  by  the  “degree  of  temporal  and/or  spatial  coherency  of  the 
wave.”  Additionally,  the  constraint  of  the  magnitudes  between  the  two  planes  results 
in  a  non-convex  set  in  the  image  space  according  to  Combettes  [6].  The  algorithm  is 
outlined  in  Figure  2.5. 
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Start 


End 


Figure  2.5:  Gerchberg-Saxton  Algorithm  [13]. 


The  inputs  to  the  algorithm  are  the  two  input  amplitude  functions  (Ref  amp) 
from  the  image  and  diffraction  planes.  A  uniformly  random  phase  generator  is  used 
to  give  the  algorithm  an  initial  starting  point.  Gerchberg  and  Saxton  offer  a  proof 
that  the  squared  error  between  the  results  of  the  Fourier  transforms  and  the  input 
amplitude  functions  must  decrease  or  remain  the  same  on  each  iteration.  They  also 
point  out  that  the  solution  reached  is  not  unique. 

Also,  in  1972,  Richardson  introduced  the  use  of  Bayesian-based  probability 
methods  to  restore  noisy  degraded  images  with  a  known  psf.  This  was  an  improvement 
over  the  Fourier  transform  methods  that  performed  well  in  a  noiseless  environment  but 
degraded  quickly  with  increased  noise.  The  degraded  image,  H ,  is  the  convolution  of 
the  original  image,  W,  with  the  point  spread  function,  S,  where  all  three  are  assumed 
to  be  discrete  probability  distribution  functions.  Two  assumptions  are  made: 


•  S  is  normalized  with  respect  to  amplitude, 

•  The  total  energy  in  H  is  equivalent  to  the  total  energy  in  W . 
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An  iterative  algorithm  is  developed  using  Bayes’  Theorem. 


The  initial  estimate  assumes  a  uniform  distribution  and  results  in  the  initial 
equation 


e  / 


Hm ,  n  Sm — i + 1 ,  n — j + 1 


_ j  /  -v  v^o  q 

%=i  n=j  2^p=a  2^q=c  ^m-p+l,n-q+l 

and  subsequent  iterations  have  the  form 


(2.11) 


where 


e  / 

\  \^ 


Wi,j,r+ 1  ^i,j,r  2_j  7  , 


Hm,nSm—i+ 1  ,n—j+\ 


m=i  n=j 


^  yb  yd  it;  .  _  sm~ 

,=i  n=i  ^-^p=a  /-^jq—c  p,q,r^m 


>,q,r^m—p-\-l,n—q+l 


a  =  (1,  m  -  K  +  l)mox; 

C  (1)  ^  L  -\~  l)?7jaa;j 

e  =  i  +  K  —  1; 


b  {jfl-i 
d  (n,  T)mjn, 

f=j  +  L-  1; 


*  =  {M};  j  =  {i,7}- 


K,  L  are  the  dimensions  of  S1^. 
I,  J  are  the  dimensions  of  Wt,r 


(2.12) 


After  the  initial  estimate,  the  W  term  on  the  right  side  of  Equation  2.12  acts  as  a  cor¬ 
rective  factor  to  the  convergence  of  the  algorithm.  The  approach  seems  to  keep  S  fixed 
throughout  the  execution  of  the  algorithm.  The  algorithm  demonstrated  convergence 
in  all  cases  tried  by  Richardson  but  no  proof  of  convergence  was  developed. 

In  1974,  Lucy  [25]  introduced  another  iterative  technique  based  on  increasing 
the  likelihood  of  the  observed  sample  each  iteration.  It  is  derived  from  Bayes’  Theo¬ 
rem  and  conserves  the  normalization  and  non- negativity  constraints  of  astronomical 
images.  The  combination  of  both  the  Richardson  algorithm  and  the  Lucy  algorithm 
is  commonly  referred  to  as  the  Richardson- Lucy  (RL)  algorithm. 

In  1982,  Fienup  [8]  compares  iterative  phase  retrieval  algorithms  to  steepest- 
descent  gradient  search  techniques  along  with  other  algorithms.  He  discusses  conver¬ 
gence  rates  and  details  some  error  functions  for  the  algorithms.  The  first  algorithm 
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discussed  is  the  error-reduction  algorithm  which  is  a  generalized  Gerchberg-Saxton 
(GS)  algorithm.  It  is  applicable  to  the  recovery  of  phase  from  two  intensity  measure¬ 
ments,  like  the  GS  algorithm,  as  well  as  one  intensity  measurement  with  a  known 
constraint.  In  the  area  of  astronomy,  the  known  constraint  is  the  non- negativity  of 
the  source  object. 

The  error-reduction  algorithm  consists  of  the  following  four  steps: 

1.  Fourier  transform  gk(x)  producing  Gk{u) 

2.  Make  minimum  changes  in  Gk{u)  which  allow  it  to  satisfy  the  Fourier  domain 
constraints  producing  G'k(u ) 

3.  Inverse  Fourier  transform  G'k(u )  producing  gk(x) 

4.  Make  minimum  changes  to  g'k  which  allow  it  to  satisfy  the  object  domain  con¬ 
straints  yielding  gk+\(x) 

where  gk(x)  and  Gk(u )  are  estimates  of  f(x)  and  F(u),  respectively. 

For  problems  in  astronomy  with  a  single  intensity  image,  the  first  three  steps 
are  identical  to  the  first  three  steps  of  the  GS  algorithm.  The  fourth  step  applies  the 
non-negativity  constraint  such  that 


9k+i(x) 


{. 9k(x )  x  £  7 k, 
0  x  E 


(2.13) 


where  7  is  the  set  of  points  dependant  on  k,  where  gk(x)  violates  the  non- negativity 
constraint,  i.e.  where  gk(x)  <  0. 

For  the  case  of  a  single  intensity  measurement  with  constraint,  Fienup’s  error 
measure  is 

a;G7  k 

He  found  that  the  “error-reduction  algorithm  usually  decreases  the  error  rapidly  for 
the  first  few  iterations  but  much  more  slowly  for  later  iterations.”  He  states  that 
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for  the  single  image  case  convergence  is  “painfully  slow”  and  unsatisfactory  for  that 
application. 

According  to  Fienup  [8]  a  method  proven  to  converge  faster  for  both  cases  is 
the  input-output  (10)  algorithm.  The  10  algorithm  differs  from  the  error-reduction 
algorithm  only  in  the  object  domain.  The  first  three  operations  are  the  same  for  both 
algorithms.  If  those  three  operations  are  grouped  together  as  shown  in  Figure  2.6 
they  can  be  viewed  as  a  nonlinear  system  having  input  g  and  output  g'.  “The  input 
g(x)  no  longer  must  be  thought  of  as  the  current  best  estimate  of  the  object;  instead 
it  can  be  thought  of  as  the  driving  function  for  the  next  output,  g'(x).” 


Figure  2.6:  Input-Output  Algorithm  [8]. 


The  basic  10  algorithm  uses  the  following  step  to  get  the  next  estimate 


9k+ i(x) 


{9k(x)  x  7*., 

9k{x)  -  Pg'k(x)  x  G  7*;. 


(2.15) 


The  output-output  algorithm  is  a  modified  10  algorithm  and  uses  the  following 
step  to  get  the  next  estimate 


9k+i(x) 


{. 9k(x )  x  £  7*, 

9k(x)  ~  P9k(x)  x  ^  Ik- 


(2.16) 
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The  hybrid  10  algorithm  uses  the  following  step  to  get  the  next  estimate 


9k+i(x) 


{. 9k(x )  x  &  7*, 

9k(x)  -  (3g'k{x)  X  G  7fc. 


(2.17) 


“The  hybrid  input-output  algorithm  is  an  attempt  to  avoid  a  stagnation  problem 
that  tends  to  occur  with  the  output-output  algorithm.”  The  output-output  algorithm 
can  stagnate  with  no  means  to  get  out  of  this  state,  but  the  hybrid  input  continues 
to  grow  until  the  output  becomes  non-negative. 

Fienup  analyzed  the  performance  differences  between  the  gradient  search  and 
the  various  10  algorithms  for  the  case  of  phase  retrieval  from  a  single  intensity  image. 
Various  values  of  (5  are  used  in  the  experiments.  The  optimal  (3  value  varied  with  the 
input  and  the  algorithm.  The  hybrid  seemed  to  have  the  best  overall  performance  but 
was  unstable  and  tends  to  move  to  a  worse  error  before  decreasing  to  a  lower  error 
than  the  other  algorithms. 

In  1988,  Ayers  and  Dainty  [1]  introduced  another  iterative  method  for  the  blind 
deconvolution  of  two  convolved  functions 


c(x) 


}{x!  )g(x 


xi)dxi . 


(2.18) 


It  is  similar  to  the  iterative  error-reduction  algorithm  but  requires  a  priori  knowledge 
of  both  functions  convolved  to  produce  the  single  convolved  image.  A  flow  chart  of 
the  algorithm  is  shown  in  Figure  2.7. 

A  priori  information  is  required  for  this  algorithm.  The  example  given  assumes 
non- negativity  of  both  functions  f{x)  and  g{x).  A  non-negative  initial  estimate  f0(x) 
is  generated  as  the  starting  point.  This  is  Fourier  transformed  to  produce  F(u)  in  Step 
1.  Then  F(u)  is  inverted  to  form  an  inverse  filter  and  multiplied  by  C(u)  to  produce 
an  estimate  G(u).  Step  3  inverse  Fourier  transforms  G(u)  to  produce  g(x).  The  image 
constraints  are  imposed  by  setting  all  points  in  g(x)  that  are  negative  to  zero  yielding 
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Figure  2.7:  Ayers  and  Dainty  Algorithm  [1], 


g(x).  In  Step  5  a  Fourier  transforms  of  g(x)  produces  G(u).  Step  6  inverts  G(u)  to 
produce  another  inverse  filter  and  multiplies  by  C(u)  to  yield  the  spectrum  estimate 
F(u).  Step  7  inverse  Fourier  transforms  F(u)  to  produce  the  function  estimate  f(x) 
and  the  non-negativity  constraint  is  applied  to  complete  the  loop. 

Two  major  problems  with  this  algorithm  are: 


1.  The  inverse  filters  in  Steps  2  and  6  are  problematic  due  to  inverting  regions  of 
low  values. 

2.  Zeros  at  particular  spatial  frequencies  in  either  F(u)  or  G(u)  result  in  voids  at 
those  frequencies  in  the  convolution. 

The  non-negativity  constraint  is  imposed  similarly  to  previous  algorithms  with 
an  additional  energy  conservation  that  seems  to  increase  the  rate  of  convergence  and 
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is  defined  as: 


fi(x) 

E 
fi(x ) 


(  fi(x),  fi(x )  >  0, 

I  0,  otherwise, 


[. fi(x)  -  fi(x)]dx, 


fi(x)  +  E/N. 


(2.19) 

(2.20) 
(2.21) 


where  N  is  the  number  of  pixels  in  the  image. 

The  Fourier  domain  constraint  is  more  complicated,  partly  to  mitigate  the  two 
significant  problems  pointed  out  previously,  and  are  summarized  as  follows: 

if  |C(«)|  <  noise  level,  then  Fi+i(u)  =  Fi(u),  (2.22) 

if  |Gi(W)|  >  \C(u)\r  then  Fi+1(u)  =  (1  -  P)Fi(u)  +  (3^-  (2.23) 

Gi(u) 

if  |Gj(u)|  <  |C(«)|,  then  — <2’24> 

Fi+1{u)  Fi(u )  C[u) 

where  0  <  (3  <  1.  On  each  iteration  the  two  Fourier  domain  estimates  are  averaged. 
The  averaging  does  not  affect  the  convergence  rate,  however,  the  convergence  rate  is 
dependent  on  j3,  but  no  method  of  finding  the  optimal  value  of  (3  was  found. 

An  additional  method  to  address  extended  regions  of  low  or  zero  values  is  in¬ 
troduced.  A  weighting  function  is  introduced  that  is  non-zero  up  to  some  band-limit. 
The  one  used  by  Ayers  and  Dainty  is  the  one  that  naturally  occurs  with  the  incoherent 
transfer  function  of  a  circular  aperture.  After  a  new  spectrum  estimate  is  produced 
by  the  averaging  process  it  is  multiplied  by  the  weighting  function.  Then,  “after  sub¬ 
sequent  imposition  of  the  image  domain  constraints  the  updated  spectrum  estimate 
is  divided  by  the  weight  function.” 

Ayers  and  Dainty  stress  “that  the  uniqueness  and  convergence  properties  of  the 
deconvolution  algorithm  are  uncertain  and  that  the  effect  of  various  amounts  of  noise 
existing  in  the  convolution  data”  is  unknown. 
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In  1990,  Seldin  and  Fienup  [39]  applied  the  Ayers  and  Dainty  (AD)  algorithm  to 
the  specific  case  of  phase  retrieval.  The  algorithm  is  compared  to  the  error-reduction 
and  hybrid  input-output  algorithms  previously  discussed  with  various  levels  of  addi¬ 
tive  Gaussian  noise.  Instead  of  implementing  the  Fourier  domain  constraint  in  Equa¬ 
tion  2.22  outlined  by  Ayers  and  Dainty,  Seldin  and  Fienup  use  a  Wiener-Helstrom 
filter  instead.  The  end  result  is  an  algorithm  that  differs  from  the  error-reduction 
algorithm  only  in  the  Fourier-domain  constraint.  The  Fourier- domain  constraint  for 
the  error-reduction  algorithm  is 

Fk(u)  =  w-SA  (2.25) 

\Fk{u)  | 

and  the  Fourier-domain  constraint  for  the  AD  algorithm  is 

Fk(u)  =  fyu)i£ML.  (2.26) 

\F„(u)  |2 

Seldin  and  Fienup  show  that  the  error-reduction  algorithm  and  the  AD  algorithm 
are  similar  both  in  form  and  in  performance.  The  error-reduction  algorithm  seemed 
to  perform  better  at  higher  noise  levels  than  the  AD  algorithm  but  the  hybrid  10 
algorithm  clearly  outperformed  both  of  them. 

In  1992,  Holmes  [18]  introduced  an  iterative  expectation-maximization  (EM) 
approach  to  the  blind  deconvolution  problem.  Primary  focus  is  on  wide-held  and 
confocal  fluorescence  microscopy.  Holmes  extends  the  maximum-likelihood  (ML)  es¬ 
timation  introduced  by  Shepp  and  Vardi  [40]  to  incoherent  imaging.  Two  constraints 
are  imposed  on  the  algorithm: 

1.  Symmetry  Constraint  -  A  circular  symmetry  based  on  the  impulse  response 
function  is  imposed  but  can  be  generalized  to  incorporate  other  constraints. 

2.  Band-limit  Constraint  -  The  knowledge  of  the  band-limit  of  the  OTF  of  the 
optical  system  is  used  to  constrain  the  solution  of  the  impulse  response. 
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Their  study  is  a  “first  attempt  at  this  type  of  approach”  with  initial  results  showing 
the  potential  use  of  ML  estimation  for  blind  deconvolution.  Extending  the  results  to 
use  MAP  estimation  or  a  penalized  log-likelihood  approach  is  recommended. 

Also  in  1992,  Ivanov  et  al.  [20]  introduced  gradient  procedures  for  the  Gerchberg- 
Saxton  algorithm.  Numerical  simulation  resulted  in: 

•  Improved  convergence; 

•  Wider  range  of  reconstruction; 

•  Better  stability  in  the  presence  of  additive  noise. 

Their  experimental  data  confirmed  the  numerical  simulation.  The  use  of  Zernike 
polynomials  was  also  introduced  for  representation  of  the  phase  function. 

Fish  et  al  [10]  developed  a  blind-deconvolution  algorithm  based  on  the  RL  algo¬ 
rithm.  The  algorithm  starts  with  an  initial  guess  for  both  the  psf  and  the  object.  The 
algorithm  then  switches  back  and  forth  between  m  iterations  to  estimate  the  next  psf 
and  m  iterations  to  estimate  the  next  object.  Better  performance  is  obtained  when  a 
priori  information  is  used  to  influence  the  initial  guess  of  the  psf. 

In  2002,  Jefferies  et  al.  [21]  discuss  the  use  of  a  constraint  on  blind  deconvolution. 
The  constraint,  in  this  case,  is  obtained  from  a  bispectrum-based  speckle  imaging 
algorithm.  The  Fourier  phase  of  the  speckle  image  is  used  to  constrain  the  recovered 
phase  estimated  by  the  blind  deconvolution  algorithm.  The  reconstructed  object 
intensity  distribution  is  modeled  as 

f(x )  =  [9(x)  (8)cu(a;)]2  (2.27) 

where  6(x)  are  the  variables  to  be  determined,  uj(x)  is  a  regularizing  function  that 
constrains  the  spatial  frequencies  of  the  reconstructed  image,  x  represents  individual 
pixels,  and  is  the  convolution  operation.  By  estimating  the  individual  9(x)  values, 
positivity  of  f(x)  is  naturally  enforced.  The  individual  pixel  psf  model  in  the  image 
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domain  is 


h(x)  =  [<y o2{x)  <3>r)(x)]  {Y>x[p2{x)  ®ry(x)]}  1  (2.28) 

where  tp(x)  are  the  variables  to  be  determined,  rj(x)  is  the  band-limited  function 

ri(x)  =  \N-2ZuP(u)ei2^u'x)/N\2  ,  (2.29) 

and  P(u)  is  the  pupil  function.  The  MAP  estimates  of  the  object  and  psf  are  obtained 
by  minimizing  the  error  developed  through  the  use  of  a  conjugate  gradient  routine. 
The  added  constraint  from  the  speckle  imaging  improves  the  convergence  properties 
of  the  blind  deconvolution  algorithm  by  narrowing  the  search  space. 

In  2004,  Likas  and  Galatsanos  [24]  introduced  a  “variational  approximation 
approach  for  Bayesian  blind  image  deconvolution.”  This  approach  can  be  viewed  as 
a  generalization  of  the  EM  algorithm.  In  an  improvement  over  an  EM  algorithm,  this 
approach  allows  the  incorporation  of  priors  for  both  the  image  and  the  psf.  It  “was  first 
introduced  in  the  machine  learning  community  to  solve  Bayesian  inference  problems 
with  complex  probabilistic  models.”  They  demonstrate  improved  performance  over 
previous  methods  through  numerical  experiments.  The  main  shortcoming  claimed  of 
the  methodology  is  that  there  is  no  analytical  way  to  evaluate  the  tightness  of  the 
variational  bound. 

2.2.2  Wave- front  sensing.  Several  other  methods  exist,  specifically  for 
detecting  the  aberrations  present  in  an  electro-optical  (EO)  imaging  system.  One 
method  is  the  Hartmann  wavefront  sensor  [29].  The  Hartmann  wavefront  sensor,  de¬ 
veloped  in  the  late  1960s,  requires  a  beam  splitter  to  divert  some  of  the  collimated 
light  in  an  imaging  system  to  a  lens  array.  Each  lens  in  the  array  produces  an  in¬ 
dividual  spot  image.  The  displacement  in  the  location  of  the  spot  image  is  used  to 
calculate  the  tilt  of  its  wavefront.  The  tilts  are  then  integrated  across  the  whole  aper¬ 
ture  to  form  the  complete  wavefront.  A  second  method  is  the  use  of  phase  diversity 
imaging  [28] .  It  is  a  technique  for  detecting  aberrations  using  a  second  diversity  im- 
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age  with  a  known  aberration.  The  known  aberration  is  usually  a  highly  calibrated 
focus  error  placed  in  the  optical  path  after  a  beam-splitter.  The  resulting  image  is 
the  diversity  image.  The  difference  between  the  conventional  image  and  the  diversity 
image  is  used  to  estimate  the  wavefront. 

A  more  recent  method  introduced  in  1988  is  curvature  sensing  [31].  Curvature 
sensing  involves  the  use  of  two  imaging  planes  equidistant  from  the  primary  imaging 
plane.  The  difference  in  irradiance  between  the  two  images  is  used  to  calculate  the 
curvature  in  the  wavefront.  The  sensitivity  of  the  curvature  sensor  in  comparison  to  a 
Shack-Hartmann  sensor  is  equivalent  [32],  This  method  requires  complicated  optics  in 
order  to  capture  both  imaging  planes  simultaneously.  For  near-simultaneous  images,  a 
membranous  mirror  can  be  used  to  change  the  focus  between  the  two  required  planes. 
In  separate  developments,  more  recent  research  resulted  in  curvature  sensing  from  a 
single  defocused  image  [17]  and  the  use  of  a  beam  splitting  cube  to  produce  the  two 
simultaneous  images  on  a  single  charge-coupled  device  (CCD)  [3]. 

All  of  these  methods  require  light  to  be  diverted  from  the  main  image  to  the 
sensor.  Diverting  light  from  the  main  image  results  in  a  lower  signal-to- noise  ratio 
(SNR)  for  the  primary  channel.  In  addition,  complicated  or  highly  calibrated  equip¬ 
ment  is  required.  An  approach  used  to  mitigate  this  problem  is  to  use  multiple  image 
frames. 

2.2.3  Multi-frame.  In  1993,  Schulz  [36]  used  ML  estimation  techniques  in 
a  multiframe  blind  deconvolution  (MFBD)  algorithm.  He  describes  an  approach  to 
blind  deconvolution  using  a  sequence  of  short-exposure  images.  Short-exposure  images 
retain  more  diffraction-limited  information  than  long-exposure  images.  However,  the 
trade  is  that  the  images  are  at  reduced  light  levels.  To  overcome  this,  a  sequence  of 
short-exposure  images  is  used  to  increase  performance.  The  intensity  image  of  each 
frame  is 

ik(y,  o,  hk)  =  ^2  hk(y\x)o(x),  y  ey  (2.30) 
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where  x  and  y  are  the  two-dimensional  coordinates  in  the  object  and  image  planes, 
respectively.  The  psf,  hk (),  is  different  for  each  frame,  k,  but  the  object  is  assumed 
to  be  the  same  for  all  frames.  Using  this  relationship  and  the  Poisson  noise  model, 
the  log-likelihood  developed  is 


K 

L(o,h)  =  Y. 

k= 1 


y ^ik{y\ o,  hk)  +  dk(y)ln(ik(y;  o,  hk)) 

y&y  y&y 


(2.31) 


Schulz  also  implements  a  penalized  ML  technique  to  avoid  convergence  to  the  trivial 
solution  of  a  point  source  as  the  estimated  image  and  the  estimates  of  the  psf  is  the 
data. 


2.2.4  Multi-channel.  No  articles  referencing  blind  deconvolution  using  mul¬ 
tiple  channels  are  currently  published.  However,  in  2004,  Cain  [5]  discusses  multichan¬ 
nel  parameter  estimation  where  two  distinct  channels  have  a  parameter  in  common. 
The  two  distinct  channels  on  the  Geostationary  Operational  Environmental  Satellite 
(GOES)  are  used  to  estimate  the  temperature  of  the  cloud  tops  that  both  channels 
observe.  The  more  accurate  temperature  estimate  leads  to  achievement  of  a  higher 
spatial  resolution  by  the  GOES.  A  similar  approach  could  be  taken  to  estimate  the 
angle  of  the  reflecting  surface  producing  the  polarization  observed  by  a  polarimeter. 
A  more  accurate  estimation  of  the  polarization  might  lead  to  a  similar  increase  in  the 
spatial  resolution  achieved  by  the  polarimeter. 

2.3  Acceleration  Techniques 

Iterative  blind  deconvolution  algorithms  are  slow.  Several  techniques  previously 
discussed  were  developed  to  accelerate  the  rate  of  convergence  of  the  iterative  algo¬ 
rithms.  However,  only  a  few  papers  were  found  that  reflect  an  attempt  to  make 
use  of  computational  techniques  such  as  parallel  implementations  or  use  of  genetic 
algorithms  to  accelerate  convergence  to  a  solution.  These  papers  are  discussed  briefly. 
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2.3.1  Parallel  Implementation.  In  [38]  Schulze  presents  parallelization  tech¬ 
niques  for  both  bispectrum  and  blind  deconvolution  and  compares  them  with  their 
serial  counterparts.  The  parallel  blind  deconvolution  implemented  is  for  multiframe 
blind  deconvolution.  The  parallelization  breaks  up  the  sequence  of  individual  frames 
and  distributes  the  individual  frames  among  additional  processors.  Using  16  frames 
of  data,  a  single  serial  image  recovery  required  87  seconds  versus  24  seconds  for  nine 
processors. 

Ingleby  and  McGaughey  [19]  present  an  algorithm  for  parallel  multiframe  blind 
deconvolution  using  wavelength  diversity.  They  use  a  sequence  of  images  that  are 
simultaneously  collected  at  M  multiple  wavelengths.  The  assumption  of  common 
path-length  errors  between  the  channels  is  exploited  to  improve  the  multiframe  blind 
deconvolution  algorithm  performance.  The  parallelization  comes  from  the  simultane¬ 
ous  collection  of  multiple  imaging  channels.  Each  individual  channel  is  processed  by 
the  multiframe  blind  deconvolution  algorithm. 

2.3.2  Genetic  Algorithms.  Schmalz  [34,35]  published  the  only  two  papers 
found  that  discuss  an  evolutionary  algorithm  (EA)  approach  to  blind  deconvolution. 
The  first  paper  is  an  initial  exploration  into  the  feasibility  of  an  EA  approach  to  blind 
deconvolution.  The  second  paper  expands  the  analysis  of  the  previous  paper’s  genetic 
algorithm  settings  versus  convergence  and  adds  an  additional  constraint  on  the  psf. 
The  performance  of  the  algorithm  is  analyzed  as  the  number  of  psf  side  lobes  expected 
is  varied. 

Both  papers  have  a  fundamental  problem  with  the  objective  function  used. 
From  [2],  “the  evaluation  process  of  individuals  in  an  evolutionary  algorithm  begins 
with  the  user-defined  objective  function, 

f:Ax^R  (2.32) 
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where  Ax  is  the  object  variable  space.”  According  to  [27],  the  objective  function  is 
“a  mathematical  statement  of  the  task  to  be  achieved.”  In  the  Schmalz  paper,  the 
objective  function  is 

g(a,a)  (2.33) 

where  a  is  the  estimate,  and  a  is  original  image.  In  the  real  world  this  is  not  possible 
since  one  does  not  have  the  original  image  to  compare  with  the  estimate.  If  one  did 
have  access  to  the  original  image,  deconvolution  is  no  longer  required. 

No  other  work  was  found  in  this  area.  EA  approaches  were  considered  during 
the  development  of  this  research  but  a  realistic  objective  function  was  not  found  to 
allow  a  realistic  solution. 

2.4  Cramer-Rao  Lower  Bounds 

Several  papers  address  various  aspects  of  CRLBs  on  atmospheric  turbulence 
parameters  and  polarimeters  related  to  this  research. 

Schulz  et  al.  [37]  discuss  CRLBs  for  estimation  of  Zcrnike  parameters  of  turbu¬ 
lence  induced  wavefront  abberations  from  conventional  and  Hartmann  sensor  images. 
The  performance  limit  for  the  estimation  of  the  first  30  Zernike  coefficients,  a,  is 
calculated  for  both  a  conventional  image  and  the  combination  of  images  from  a  Hart¬ 
mann  sensor  array.  The  collected  data,  T> ,  for  all  the  detector  elements  is  the  set 
T>  =  {d[p]}  where  p  is  the  two-dimensional  coordinates  of  the  data/image  and  the 
conditional  probability  mass  function  of  the  data  is 


P(V-a )  =  Y[e{-h[p'a]\h[p]a])dlp]/d[p}\, 

p 

where  the  psf  for  the  individual  images  is 


J  wp(y)h(y,  a)dy. 


(2.34) 


h(pi  a) 


(2.35) 


wp(y )  is  a  “normalized  weighting  function  that  models  the  spatial  region  of  interest 
for  the  pth  detector.”  The  individual  elements  of  the  Fisher  information  matrix  are 


[./(o)  nm 


d2  In  P(T>]  a) 

0(ln()(lrn 


(2.36) 


and  the  Cramer-Rao  lower  bound  for  any  unbiased  estimator  of  an  unknown  Zernike 
parameter  vector  a  is 

E[(a  -a)(a-  a)T]  >  [J(a)]"1  (2.37) 

where  the  inequality  is  understood  to  be  entry-wise.  However,  since  a  is  a  random 
parameter  vector  the  resulting  lower  bound  is  also  random.  Their  results  are  the 
average  of  200  realizations  and  the  bound  values  are  normalized  by  the  expected 
number  of  photons. 

Gamiz  and  Belsher  [11, 12]  discuss  the  effects  of  detection  noise  and  calibra¬ 
tion  errors  on  the  performance  of  a  four-channel  polarimeter.  Their  work  “provides 
a  framework  for  developing  polarimeter  designs  with  minimal  noise  sensitivity  and 
evaluating  the  noise  performance  of  existing  designs.” 

This  chapter  presented  several  topics  from  the  description  of  a  basic  polarime¬ 
ter  to  the  development  of  approaches  to  blind  deconvolution  using  standard  electro- 
optical  imaging  systems.  Additionally,  previous  CRLB  calculations  related  to  Zernike 
coefficient  estimation  and  the  effects  of  noise  and  calibration  errors  on  polarimeters 
were  presented.  The  development  of  the  spatial  resolution  CRLB  of  a  basic  two- 
channel  polarimeter  is  presented  in  the  next  chapter. 


29 


III.  Polarimeter  Cramer-Rao  Lower  Bounds  on  Spatial 

Resolution 


This  research  utilizes  statistical  models  to  develop  Cramer-Rao  lower  bounds 
(CRLB)  for  the  resolution  of  a  polarimetric  sensor  and  a  non-polarization  sen¬ 
sitive  imaging  system.  The  first  of  the  two  cases  consists  of  a  single  image  of  two 
point  sources  and  is  used  for  comparison  with  the  second  case.  The  second  case 
models  two  simultaneous  images  that  would  result  from  the  incoming  beam  of  light 
passing  through  a  polarizing  beam  splitter.  Both  cases  assume  an  ideal  unaberrated 
image  collection  with  no  atmospheric  turbulence.  The  Fisher  information  (FI)  matri¬ 
ces  are  then  formed  from  the  individual  elements.  The  CRLBs,  which  are  the  lower 
bounds  on  the  variances  for  the  parameters  estimated,  are  then  produced  from  the 
FI  matrix.  The  CRLBs  are  developed  in  three  steps.  The  first  is  the  development 
of  the  spatial  resolution  CRLB  of  a  standard  imaging  system  with  a  single  channel 
without  turbulence.  The  second  step  is  to  extend  the  spatial  resolution  CRLB  to  a 
two-channel  polarimeter.  The  final  step  is  to  introduce  an  atmospheric  turbulence 
model  based  on  a  Zcrnike  polynomial  phase  screen.  This  is  followed  by  a  comparison 
of  the  polarimeter  spatial  resolution  CRLB  to  that  of  a  standard  imaging  system. 


3.1  Bound  for  Unpolarized  Primary  Channel 

The  base  case  for  this  analysis  uses  two  point  sources  as  the  image.  The  pa¬ 
rameters  estimated  are  the  amplitudes  of  the  two  point  sources  (o o  and  oi)  and  the 
separation  distance  (A).  The  CRLBs  are  then  produced  from  the  FI  matrix  developed 
from  the  image  model. 


3.1.1  Image  model.  The  following  probability  mass  function  (PMF)  is  used 
for  development  of  the  single  image  used  in  this  scenario 


PMF(  a,O0,Oi) 


nn 

x  y 


i(x,  y'\d(x’y') e~dx'y) 
d(x,  y)\ 


(3.1) 
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where  d(x,y )  is  the  data  or  image  captured  and  i(x,y )  is  the  noiseless  image  defined 


as 

i(x,v)  =  ^2^2h(x- z,y-w)o(z,w).  (3.2) 

z  w 

The  term  o(z,  w)  is  defined  as 

o(z,  w)  =  0q5(z,  w)  +  o\S(z  —  A,  w),  (3.3) 


where  5  is  the  Kronecker  delta  function,  and  Oo  and  o\  are  the  intensities  of  the 
two  point  sources  in  the  image.  The  value  A  is  the  separation  between  the  points. 
Therefore, 


i(x,y)  =  EE  h(x  —  z,y  —  w)(oqS(z,  w )  +  o\5{z  —  A,  w)). 


(3.4) 


The  point  spread  function  [15]  psf  chosen  for  this  analysis  is  a  radial  symmetric 
Gaussian  centered  at  the  origin  (0,  0)  with  variance  a2  >  0, 


h(x,y)  = 


27 ra2 


exp 


■( x 2  +  y2 
2  u2 


(3.5) 


It  is  chosen  for  ease  of  computation  of  its  derivative  when  shifted  with  respect  to  A. 
The  derivative  of  the  shifted  psf  is 


r  //  an  <9  ,  ,  .  .  —x  +  A 

h  (x  -  A,  y)  =  —  h(x  -  A,  y)  =  exp 


-((x-A)2  +  1/2) 
2cr2 


(3.6) 


and  is  required  to  compute  certain  derivatives  of  the  image  with  respect  to  A. 
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The  log-likelihood  function  L (L0,  Li,  L2)  is  the  log  of  the  PMF.  For  the  three 
parameter  case  in  this  scenario  L0  =  A,  Li  =  oq,  and  L2  —  0\. 


L(A,  Oq,  O i) 


Kim 


x,y)„-i(x,y) 


x  y 


i(x,  y)d^'y,e 
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^2^2[d(x,y)  In (i(x,y)) 

x  y 

^2^2[d(x,y)  In (i(x,y)) 


x  y 


) 

i(x,y)  -  ln(d(x,y)!)] 

i{x,y)\  -  5Z5Zln(d(®,J/)!)-  (3-7) 

x  y 


3.1.2  Fisher  Information  Matrix  and  CRLB.  The  CRLB  for  the  individual 
coefficients  estimated  are  [45] 


CT; 


>  r  =  [j 


-n 


(3.8) 


where  J  is  the  Fisher  information  matrix  and  is  defined  as  follows: 
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(JL.dL, 


(3.9) 


(3.10) 


where  E  is  the  expected  value  and 
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(3.11) 
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—E 
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(3.12) 


Therefore 


J«  =  EE 


a;  y 


1  di(x,y)  di(x,y) 
i(x,y)  dLi  dLj 


(3.13) 


3.1.3  Derivatives  of  the  Image.  In  order  to  calculate  the  derivatives  of  the 
image,  one  additional  definition  [7]  is  required.  Let  /  be  a  continuous  function  defined 
on  the  interval  Let  tQ  £  [ti,f2]-  Then 


f't2 


f(t)S’(t  -  tc)dt  =  (— l)/'(to). 


(3.14) 


The  Fisher  information  matrix  is  then  populated  using  the  following  three  image 
derivatives: 


di(x'y)  4ee  h{x  —  z,y  —  w)(o0S(z,  w)  +  0\5{z  —  A,  w)) 


dA 
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(3.15) 
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Z  W 

EE  h(x-  z,y-w)8(z, 

Z  W 

h(x,y )  (3.16) 

d 

£-EE  /i(x  —  z,y  —  w)(o05(z ,  it;)  +  0i5(z  —  A,  it;)) 

^  2: 

E  E  h(x  ~ z> y  ~  w^(z  -  A’ ^ 

2:  w 

h(x-A,y )  (3.17) 


Using  the  previous  results  and  Equation  3.14  the  elements  of  the  Fisher  information 
matrix  are 
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(3.18) 
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(3.19) 
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(3.21) 


(3.22) 


(3.23) 

(3.24) 

(3.25) 


(3.26) 


3.1.4  Cramer-Rao  Lower  Bound.  The  CRLBs  or  smallest  possible  variances 
of  the  estimated  parameters  are  calculated  by  taking  the  inverse  of  the  FI  matrix.  The 
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individual  bounds  are  located  along  the  diagonal  of  the  resulting  matrix. 


Jl°  =  J~l  = 


<4 

J  A°o 
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JO0°l 

JOlA 

JOlOQ 

4 

Alternatively,  the  individual  bounds  can  be  found  as  follows 


<4  =  det 
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OlOl 


det[J] 


ao1  = 


J A  A  J Aoo 

JoqA 
0 
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det[J] 


(3.27) 


(3.28) 


(3.29) 


(3.30) 


3.2  Polarimeter  Bounds  without  Atmospheric  Distortion 

Similar  to  the  base  case,  two  point  sources  comprise  the  image  modeled.  The 
parameters  are  expanded  to  include  intensities  of  the  point  sources  to  be  estimated 
in  two  polarized  images.  The  scenario  assumes  a  perfect  PBS  to  produce  the  two 
images.  Thus  the  parameters  estimated  are  the  amplitudes  of  the  two  point  sources 
( Ooh,Oov,Oih  and  Oi„)and  the  separation  distance  (A).  The  CRLBs  are  then  produced 
from  the  FI  matrix  developed  from  the  image  model  in  a  similar  manner  to  those  in 
the  base  case.  Additional  discussion  looks  at  the  impact  of  a  less  than  perfect  PBS. 


3.2.1  Image  Model.  The  following  PMF  is  used  for  development  of  the 
CRLB  for  a  scenario  using  a  sensor  that  captures  two  simultaneous  images  using  the 


36 


two  different  polarity  channels  produced  by  a  PBS: 


P M F^  O0h  O0v  ,0\h,01v)  |  | 


ih(x,  y')dh(x’y)e~ih(x,v^  iv(x ,  yYvP,v^ e~ivP'y^ 


x  y 


dh(x,y)\ 


dv(x,  y)\ 


(3.31) 


The  image  dehnitions  for  notation  purposes  arbitrarily  define  one  image  as  the  hori¬ 
zontally  (h)  polarized  image  and  the  other  as  the  vertically  (v)  polarized  image.  The 
individual  images  are  thus  defined,  similarly  to  the  base  case,  as 


ih(x,y)  =  EE  h(x  -  z,y  -  w)(ooh5(z,w )  +  olh5(z  -  A ,w))  (3.32) 

Z  W 

iv{x,y )  =  ^h(x-  z,y-w)(o0v5(z,w)  +  olv8(z  -  A ,w))  (3.33) 

2  W 

The  log-likelihood  function  for  the  dual  polarity  case  is  L(A,  o^h,  Oqv,  o\h,  Oiv).  By 
development  similar  to  that  used  for  the  base  case,  the  result  is: 


L(A,  o0h,  o0v,  oih,  olv)  =  In (PMF) 

«EEirf  h(x,  y )  In (ih(x,  y ))  -  ih{x,  y)  +  dv(x,  y)  ln(i^(x,  y ))  -  iv(x,  y)] 

x  y 


3.2.2  Polarimeter  Fisher  Information  Matrix  and  CRLB.  The  Fisher  infor¬ 
mation  matrix  is  composed  of  the  following  elements 
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(3.34) 


(3.35) 
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3.2.3  Polarimeter  Image  Derivatives.  The  individual  derivatives  of  the 
image  with  respect  to  the  parameters  estimated  are 
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(3.44) 

3.2.4  Individual  FI  Elements  for  Polarimeter.  The  individual  elements  of 
the  FI  matrix  J  are  calculated  using  the  image  derivatives  of  the  previous  section  and 
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3.2.5  FI  Matrix  and  Degree  of  Polarization.  The  degree  of  polarization 
P  affects  the  four  intensities  used  in  the  image  model  and  the  calculations  of  the 
CRLBs.  The  value  of  P  ranges  from  0  for  completely  unpolarized  light  to  1  for 
completely  polarized  light.  With  P  <  1  the  resulting  FI  matrix  is 
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For  a  completely  polarized  image,  with  one  point  source  horizontally  polarized  and 
the  other  vertically  polarized,  four  more  terms  in  the  FI  matrix  are  zero.  This  is 
due  to  two  of  the  four  intensity  terms  going  to  zero  in  the  limit.  For  instance,  let 
Ooh  =  Oiv  —  1  and  o0v  =  oih  =  0,  then  the  resulting  FI  matrix  is 
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3.2.6  Cramer-Rao  Lower  Bound  of  Polarimeter.  Just  like  the  base  case, 
the  CRLBs  are  found  by  inverting  J. 
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(3.63) 


3.2.7  Analysis.  Due  to  the  cumbersome  mathematics  involved  with  a  direct 
comparison  between  the  two  cases,  a  Mat  lab®  analysis  is  developed  and  the  results 
used  to  evaluate  any  performance  differences.  In  order  to  evaluate  the  FI  matrix  a 
generic  telescope  with  an  aperture  size  of  one  meter  and  no  centFinally,  some  possible 
extensions  of  this  analysis  are  proposed. ral  obscuration  is  used  to  simulate  viewing 
the  two  point  sources.  For  the  polarimeter  images,  the  effect  of  a  polarizing  beam 
splitter  is  modeled.  The  spacing  in  the  detector  plane  is  1/iV  where  N  is  the  number 
of  pixels  along  the  side  of  the  square  detector  array.  For  the  analysis  in  this  paper 
N  =  256.  Another  necessity  involves  the  use  of  a  pseudo- inverse  technique  to  invert 
the  FI  matrix. 


Figure  3.1  shows  a  graphical  view  of  the  simulation  during  operation  with  a  A  of 
20  pixels  separating  the  two  point  sources.  The  point  sources  have  been  made  larger 
in  order  to  be  visible  in  the  figure,  but  during  the  actual  analysis  they  are  reduced  to 
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Figure  3.1:  Visualization  of  Running  Simulation. 

a  standard  deviation  for  the  Gaussian  of  approximately  one  pixel  width  to  simulate 
the  delta  function.  The  plots  of  the  bound  calculated  for  A  are  shown  in  Figure  3.2 
as  the  polarization  is  varied. 

For  the  particular  run  shown,  the  base  case  value  for  <7  a  is  0.07  while  the  results 
for  the  bound  with  p  =  1  is  a  a  =  0.035.  Therefore,  under  ideal  conditions,  the  bound 
for  two  oppositely  polarized  sources  is  approximately  50%  of  the  bound  for  the  base 
case. 


3.2.8  Transmission  and  Reflection  Efficiency.  A  quick  sampling  of  polar¬ 
izing  beam  splitters  showed  general  transmission  and  reflection  efficiencies  of  better 
than  95  and  99%  ,  respectively.  These  values  are  added  into  the  modeling  as  efficiency 
parameters  with  the  result  that  the  bound  increases  slightly.  Using  more  realistic 
transmission  and  reflection  efficiencies  results  in  an  approximately  4%  increase  in  the 
bound. 
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Figure  3.2:  CRLB  on  estimation  of  pixel  separation  versus  actual  pixel  separation. 


The  results  of  this  experiment  show  that  the  resolution  of  two  point  sources  is 
doubled  when  the  sources  are  polarized  90  degrees  from  each  other.  The  performance 
degrades  as  the  degree  of  polarization  decreases  from  1  to  0.  At  p  =  0  the  images 
collected  are  completely  unpolarized  and  the  resulting  bound  is  the  same  as  for  the 
base  case  of  a  single  image.  For  realistic  efficiencies  for  the  PBS,  the  resulting  bound 
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for  unpolarized  sources  is  approximately  three  percent  lower  than  the  base  case.  The 
decrease  in  the  bound  or  increase  in  resolution  is  directly  tied  to  the  amount  of  po¬ 
larization  of  the  source  imaged.  Therefore,  potential  does  exist  to  improve  resolution 
through  the  development  of  algorithms  that  can  make  use  of  polarization  data  from 
captured  images. 

3.3  Spatial  Resolution  Bound  with  Atmospheric  Turbulence 

The  following  section  extends  the  two  channel  polarimeter  case  to  include  the 
effects  of  atmospheric  distortion.  Similar  to  the  previous  development,  two  point 
sources  comprise  the  image  modeled.  The  parameters  are  expanded  to  include  the 
individual  Zernike  coefficient  parameters.  The  scenario  still  assumes  a  perfect  PBS 
to  produce  the  two  images.  Thus  the  parameters  estimated  are  the  amplitudes  of 
the  two  point  sources  (ooh,oov,oih  and  oiv),  the  separation  distance  (A),  and  a  set  of 
Zernike  parameters  (a2,a3,...,an_i,an).  The  average  CRLBs  are  then  produced  from 
the  averaging  of  multiple  instantiations  of  the  FI  matrix  developed  from  the  image 
model  and  then  inverting  it. 

3.3.1  Image  Model.  The  PMF  used  in  Equation  3.31  is  extended  to  include 
the  atmospheric  turbulence  model 

PMF(  A,ooh  ,o0„  ,0lh,01v  ,0,2 ,0.3,-  ,a„-i,a„) 

Xfi(x  y  od^bP’^e  ^  y  ) 0 

dh(x,y)\  ^  ^  dv(x,y)\ 

The  image  definitions  for  notation  purposes  arbitrarily  define  one  image  as  the  hori¬ 
zontally  (h)  polarized  image  and  the  other  as  the  vertically  (v)  polarized  image.  The 


)  (3.64) 
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individual  images  are  thus  defined,  similarly  to  the  base  case,  as 


i(x,y,a) 

=  ih(x,  y;  a)  +  iv(x,  y;  a) 

=  EE  h(x  -  z,y  -w;  a)((ooh  +  o0v)5(z,w)  +  (olh  +  oiv)S(z  -  A ,w)) 

Z  W 

(3.65) 


where  h(x,  y;  a)  is  the  psf  conditioned  on  the  turbulence. 

3.3.2  The  Point  Spread  Function.  To  account  for  atmospheric  turbulence 
and  the  effect  of  the  imaging  system  in  a  closer  approximation  to  reality  a  more 
accurate  model  of  the  psf  than  that  presented  in  Section  3.1.1  is  required.  With  the 
assumption  of  an  abberation  free  imaging  system,  the  only  contribution  to  the  psf  is 
the  pupil  itself.  When  the  point  source  illumination  is  sufficiently  far  away  the  field 
arrives  at  the  aperture  as  a  plane  wave.  Under  these  conditions  the  image  of  the 
point  source  is  the  Fourier  transform  of  the  field  at  the  aperture  masked  by  the  pupil 
function  and  multiplied  by  the  wave  abberation  function  caused  by  the  atmosphere. 
This  yields  the  psf 
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and  the  psf  for  the  shifted  point  is 
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(3.67) 


Thus  h(x  —  A ,y;a)  is  equal  to  h(x,y,a )  shifted  by  A  using  the  Fourier  transform 
shift  theorem. 


3.3.3  Log-likelihood  Function.  The  log-likelihood  function  for  the  dual  po¬ 
larity  case  is  found  by  taking  the  log  of  the  PMF  as  shown  in  Eqn.  3.64.  It  is  a 
function  of  all  the  parameters  estimated  in  the  model. 
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x  y 

+dv(x,  y )  In (iv(x,  y,a))~  iv(x,  y,  a)]  (3.68) 


3.3.4  Polarimeter  Fisher  Information  Matrix  and  CRLB.  The  Fisher  in¬ 
formation  matrix  is  composed  of  the  expected  values  of  the  partial  derivatives  of  the 
log- likelihood  function  [45].  The  individual  elements  in  the  FI  matrix  are 
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The  resulting  FI  matrix  is  shown  in  Equation  3.70. 


(3.69) 
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The  Zernike  polynomials  are  an  orthogonal  basis  set,  therefore  the  cross  terms 
are  zero  and  the  FI  matrix  simplifies  to  Equation  3.71. 
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3.3.5  Polarimeter  Image  Derivatives.  The  individual  derivatives  of  the 
image  with  respect  to  the  parameters  estimated  are 
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(3.72) 
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3.3.6  Individual  FI  Elements  for  Polarimeter.  The  individual  elements  of 
the  FI  matrix  J  are  calculated  using  the  image  derivatives  of  the  previous  section  and 
are: 
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3.3.7  Partial  Derivatives  of  the  psf.  The  following  partial  derivatives  of  the 
psf  are  required  for  calculation  of  the  individual  FI  elements. 
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The  partial  derivatives  of  the  psf  with  respect  to  the  individual  Zernike  coeffi¬ 
cients  are 
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5.5.5  CRLB.  Using  the  previous  partial  derivatives,  the  FI  matrix  is  filled 
in  and  the  inverse  is  calculated.  The  resulting  diagonal  elements,  shown  in  Equation 
3.116,  are  the  CRLBs  for  the  parameters  estimated.  The  primary  element  studied 
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in  the  following  sections  is  the  first  diagonal  element  a\.  This  is  the  CRLB  on  the 
estimate  of  the  pixel  separation  A. 
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3.3.9  Effect  of  Atmospheric  Turbulence  on  CRLB  .  Using  the  expanded 
FI  matrices,  the  average  lower  bound  on  the  estimate  of  the  pixel  separation  and 
other  parameters  are  calculated.  This  is  accomplished  by  generating  a  new  set  of 
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Zernike  coefficients  for  each  change  in  parameters.  Multiple  iterations  are  then  av¬ 
eraged  together  to  yield  the  average  CRLBs  for  the  parameters  estimated.  In  order 
to  compute  an  analytical  solution,  some  additional  parameters  are  required.  The 
following  parameters  are  used  for  generation  of  most  of  the  graphs  in  this  analysis: 

•  Aperture  size  of  1  meter 

•  Fried  parameter  (rD)  of  1  meter 

•  Focal  length  of  10  meters 

•  Object  intensity  of  O.lmW 

•  Grid  size  of  32  by  32  pixels 

•  Wavelength  of  800nm 

3.3.10  Normal  Imaging  System.  This  section  discusses  the  average  lower 
bound  on  the  estimate  of  the  pixel  separation  for  a  normal  (i.e.  non-polarimetric) 
sensor.  The  number  of  Zernike  coefficients  is  varied  and  the  impact  on  the  average 
lower  bound  is  shown  in  Figure  3.3. 


Figure  3.3:  Average  bound  relative  to  pixel  separation. 

The  primary  focus  of  this  research  is  on  the  CRLB  of  the  estimated  pixel  sepa¬ 
ration  of  the  two  point  sources.  Figure  3.3(a)  shows  the  average  CRLB  as  the  actual 
pixel  separation  in  the  model  is  changed.  The  lowest  line  is  the  basic  model  without 
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any  Zernike  coefficients  estimated.  The  upper  grouping  of  lines  are  for  estimates  when 
two,  four,  six,  and  eight  coefficients  are  added  to  the  model.  Figure  3.3(b)  shows  the 
standard  deviations  of  the  data  used  to  generate  the  bound  plots  to  the  left.  The  stan¬ 
dard  deviation  for  the  model  without  any  Zernike  coefficients  is  on  the  bottom  of  the 
plot  at  values  in  the  10~14  range.  Thus,  introducing  atmospheric  turbulence  into  the 
model  increases  the  lower  bound.  The  bound  follows  a  similarly  shaped  curve  as  the 
non-turbulent  model.  All  plots  asymptotically  approach  some  value  as  the  separation 
gets  larger.  The  determination  of  that  value  is  discussed  later  in  this  chapter. 

3.3.11  Polarimeter  Imaging  System.  Now  that  the  behavior  for  the  normal 
imaging  system  is  known,  the  behavior  of  the  polarimeter  model  is  determined.  Keep¬ 
ing  the  same  parameter  settings  as  the  normal  system,  the  polarimeter  lower  bound 
on  the  estimation  of  A  is  plotted  against  the  actual  pixel  separation  with  the  addition 
of  various  induced  polarization  states  of  the  point  sources. 


(a)  Polarimeter  without  Zernike  Coefficients.  (b)  Polarimeter  with  Zernikes  2  and  3. 


Figure  3.4:  Average  polarimeter  bounds  relative  to  pixel  separation. 


Figure  3.4(a)  shows  the  effect  of  polarization  on  the  average  lower  bound  on 
the  estimated  pixel  separation  as  the  degree  of  polarization  goes  from  0  to  close  to 
1.  A  note  on  the  model  regarding  a  completely  polarized  source  is  necessary  here. 
With  a  degree  of  polarization  of  1,  a  discontinuity  is  reached  and  data  is  erroneous. 
Therefore,  a  degree  of  polarization  of  0.99  is  chosen  to  show  the  average  bound  close 
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to  p  —  1.  One  should  also  note  that  as  the  degree  of  polarization  increases,  the  bound 
flattens  out  until  the  limit  on  the  bound  is  reached  for  estimating  the  location  of  a 
single  point  source. 

In  Figure  3.4(b),  the  estimates  of  the  first  two  Zernike  coefficients  (piston  is 
ignored)  related  to  tilt  are  introduced.  The  bounds  have  a  relative  shift  upward. 
The  top  line  in  the  plot  corresponds  to  the  unpolarized  source  case  and  follows  the 
grouping  of  lines  seen  in  Figure  3.3(a)  for  the  normal  imaging  system.  As  the  degree  of 
polarization  increases,  the  average  lower  bound  decreases.  The  limit  is  approached  as 
the  degree  of  polarization  approaches  1.  Based  on  a  definition  of  minimum  resolving 
power  for  the  system  located  at  the  point  where  its  estimate  equals  the  actual,  Figure 
3.4  shows  that  polarization  effects  can  improve  resolving  power.  For  this  particular 
set  of  parameters,  the  improvement  is  approximately  20%. 

3.3.12  Other  System  Parameters.  The  effects  of  other  parameters  on  the 
spatial  frequency  lower  bound  and  other  estimates  are  explored  in  this  section.  First, 
the  effect  the  number  of  Zernike  coefficients  estimated  have  on  the  spatial  frequency 
bound  and  on  source  intensity  estimates  are  explored.  Next,  all  parameters  are  fixed 
and  just  the  source  intensities  are  varied  to  explore  its  effect  on  the  bound.  Finally,  a 
very  brief  look  is  taken  at  the  effect  of  the  Fried  parameter  on  the  model  to  see  how 
increasing  the  turbulence  in  the  model  affects  the  results.  It  is  fixed  at  eight  times 
smaller  than  the  aperture  and  the  impact  of  partial  compensation  by  an  AO  system 
is  explored. 

An  understanding  of  how  the  model  changes  as  the  number  of  Zernike  coefficients 
estimated  increases  is  explored.  It  is  noted  that  the  average  lower  bound  on  the 
estimate  of  A  increases  as  the  number  of  coefficients  increases.  In  order  to  quantify 
the  change,  all  the  parameters  are  fixed.  The  average  pixel  separation  bound  on 
estimatation  is  then  calculated  at  an  actual  pixel  separation  of  five  pixels.  Due  to  the 
fairly  quick  execution  time,  an  arbitrarily  large  set  of  two  hundred  randomly  generated 
Zernike  coefficients  were  used  to  generate  the  plots  shown.  They  are  calculated  for 
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both  the  normal  imaging  system  and  the  polarimeter  at  each  polarization  setting  for 
every  additional  Zernike  coefficient  up  to  85. 

Figure  3.5(a)  plots  the  result  for  the  normal  imaging  system.  Figure  3.5(b) 
plots  the  results  for  the  polarimeter  for  each  of  the  polarization  settings.  As  the 
degree  of  polarization  increases  the  plot  lines  decrease  but  retain  similar  shapes.  The 
plot  for  the  polarization  line  close  to  1  gets  more  erratic  as  the  degree  of  polarization 
gets  closer  to  1.  This  is  due  to  the  discontinuity  in  the  model  when  the  degree  of 
polarization  is  1.  All  of  the  lines  are  similar  in  shape  to  the  plot  in  Figure  3.5(c). 
This  is  the  cumulative  sum  of  the  expected  values  for  the  covariance  of  the  Zernike 
coefficients. 


(a)  Normal  average  pixel  separation  bound  (b)  Polarimeter  average  pixel  separation 

bound 


Number  of  Zernike  coefficients 

(c)  Zernike  coefficient  covariance  cumulative 
sum 

Figure  3.5:  Average  bound  relative  to  number  of  Zernike  coefficients  estimated. 
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Although  the  average  spatial  frequency  CRLB  changes  as  the  number  of  Zernike 
coefficients  changes,  the  average  CRLB  of  the  point  source  intensity  estimates  change 
very  little.  Very  slight  variations  of  the  bounds  are  noted  in  Figure  3.6,  but  the 
plots  are  almost  flat  as  the  number  of  coefficients  is  increased.  When  the  model  is 
unpolarized  all  estimates  for  the  point  sources  have  the  same  values.  As  the  degree 
of  polarization  increases,  the  intensities  of  0oh  and  0\v  increase  in  the  model  with  a 
corresponding  increase  in  the  average  lower  bound  of  their  estimates.  The  converse 
is  true  for  0ov  and  Oih. 

xIO"5  xIO'5 


6 


P=0.00 


Number  of  Zernike  coefficients  estimated 


(a)  Estimates  of  Ooh  (b)  Estimates  of  Oov 


(c)  Estimates  of  Olh 


(d)  Estimates  of  Olv 


Figure  3.6:  Intensity  estimates  vs.  number  of  coefficients  estimated. 


Another  parameter  that  affects  the  model’s  estimates  is  the  light  intensity  of 
the  point  sources.  For  simplicity  both  source  are  fixed  at  the  same  magnitude.  While 
all  other  parameters  are  fixed,  the  intensity  is  adjusted  to  produce  the  curve  shown 
in  Figure  3.7.  As  the  intensity  is  increased,  the  average  spatial  frequency  CRLB 
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decreases  asymptotically  to  one  pixel  width.  Conversely,  as  the  intensity  is  decreased, 
the  average  spatial  frequency  CRLB  increases. 


Figure  3.7:  Average  bound  vs.  point  source  intensity. 

The  preceding  analysis  in  this  paper  assumes  a  D /r0  fraction  of  1,  where  ra  is 
the  Fried  parameter  (ra)  [33].  This  section,  however,  discusses  the  effect  of  an  rQ  value 
less  than  the  aperture  size  such  that  D/ra  =  8.  The  plots  shown  in  Figure  3.8  reflect 
this  increased  “turbulence.”  The  left  image,  Figure  3.8(a),  plots  the  average  spatial 
frequency  lower  bound  for  each  polarization  case.  The  results  are  from  the  averaging 
of  200  independent  iterations  at  each  polarization  setting. 

The  impact  of  an  AO  system  is  also  modeled.  Assuming  an  AO  system  can 
remove  95%  of  the  tilt  coefficients,  the  results  shown  in  Figure  3.8(b)  are  attained. 
The  initial  part  of  the  plots  is  greatly  smoothed  out  and  the  overall  bounds  are 
improved  significantly.  For  the  last  case,  the  following  characteristics  are  assumed 
regarding  the  “AO”  system’s  compensation  of  abberations: 

•  95%  tip/tilt  compensated 

•  80%  defocus  and  astigmatism  compensated 

•  50%  coma,  trifoil,  and  spherical  compensated 
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Figure  3.8:  Effect  of  removing  lower  order  abberations. 

The  bounds  are  significantly  smoothed  out,  in  Figure  3.8(c),  for  the  entire  range  of 
Zernike  coefficients.  The  results  show  that  the  model  can  incorporate  AO  parameters 
for  estimation  of  the  spatial  frequency  CRLB  under  different  atmospheric  conditions. 
Thus,  if  one  knows  the  performance  characteristics  of  an  AO  system  and  the  rQ  value 
during  the  observation,  one  can  use  the  model  to  predict  the  lower  bound  on  spatial 
frequency  resolution  that  can  be  attained  by  any  algorithm. 

This  chapter  incrementally  developed  the  spatial  resolution  CRLB  for  a  two- 
channel  polarimeter  with  atmospheric  turbulence.  Starting  with  a  normal  imaging 
system.  Then  the  PMF  is  extended  to  two  channels.  Finally,  the  psf  is  extended  to 
include  a  Zernike  polynomial  based  phase  screen  turbulence  model.  The  CRLB  on 
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the  estimate  of  the  pixel  separation  is  compared  as  the  models  increase  in  complexity. 
The  next  chapter  details  the  development  of  calibrated  test  data  used  to  test  any 
algorithms  developed,  along  with  a  new  approach  for  detecting  focus  aberrations  in  a 
single  imaging  channel. 
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IV.  Laboratory  Data  Calibration 


This  chapter  discusses  the  data  used  in  development  and  testing  of  the  polarime- 
ter  multi- frame  blind  deconvolution  (PMFBD)  algorithm.  Since  tip/tilt  aber¬ 
rations  do  not  affect  short-exposure  image  quality  [33]  and  are  fairly  easy  to  remove 
through  other  methods  they  are  not  used  for  algorithm  development  for  this  research. 
The  defocus  aberration  used  in  algorithm  development  is  one  of  the  next  most  signif¬ 
icant  aberrations  caused  by  atmospheric  turbulence.  It  is  also  a  simple  aberration  to 
create  in  the  laboratory  environment  by  simply  moving  the  lens.  This  allows  testing 
of  the  developed  algorithm  with  calibrated  data  in  addition  to  simulated  data. 

Following  the  section  on  the  focus  aberration,  an  innovative  focus  aberration 
detection  (FAD)  algorithm  is  presented.  The  algorithm  is  used  to  detect  the  magni¬ 
tude  of  the  focus  aberration  present  in  a  single  imaging  channel.  Due  to  the  multiple 
optical  paths  in  a  polarimeter,  focus  aberration  detection  is  important  in  calibration 
of  the  separate  channels.  The  results  of  the  FAD  algorithm  are  used  for  validation 
of  the  laboratory  setup  and  simulation  used  to  produce  additional  images  for  testing. 
The  chapter  concludes  with  a  brief  discussion  of  a  potential  autofocus  algorithm  using 
the  FAD  algorithm  and  possible  methods  of  improving  the  FAD  algorithm’s  wallclock 
execution  time. 


4-1  Focus  Aberration 

There  are  many  factors  that  can  contribute  to  an  out  of  focus  image.  From 
geometrical  optics,  equations  for  focus  are  well  established.  Using  simplifying  tech¬ 
niques,  a  basic  optical  system  can  be  simplified  to  an  object,  a  lens,  an  image  plane, 
and  the  distances  between  them.  Using  this  very  simple  system,  the  location  of  the 
image  plane  is  determined  by  the  equation  [16] 


^image  d0bject  •  f  /  (^d  object,  .[) 


(4.1) 


where  d0bject  is  the  distance  between  the  object  of  interest  and  the  lens,  dimage  is  the 
distance  between  the  lens  and  the  image  plane,  and  /  is  the  focal  length  of  the  lens. 
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Regardless  of  the  cause  of  the  focus  aberration,  its  characteristic  shape  remains 
the  same.  A  Zernike  polynomial  representation  is  used  to  create  the  phase  screen  in 
the  algorithm,  where  the  phase  screen  is 


phase  screen 


pupil  ■  exp 


— ja  ■  —  ■  Zernike 4 

Z7T 


(4.2) 


For  simplicity,  a  circular  pupil  is  chosen,  a  is  the  magnitude  of  the  focus  aberration 
in  waves  and  Zernike 4  [33]  is  the  Zernike  polynomial  for  focus.  By  adding  the  factor 
A / 27 r,  where  A  is  the  wavelength  of  light  in  the  scenario,  the  analysis  becomes  wave¬ 
length  independent  [4],  The  phase  screen  is  then  used  to  generate  the  OTF  used  in 
the  deconvolution  algorithm. 


The  OTFs  for  several  focus  aberrations  are  shown  in  Figures  4.1(a),  4.1(b),  and 
4.1(c).  The  aberrations  are  normalized  for  incoherent  imaging  such  that  the  corners 
of  the  plots  are  at  fx  =  2  fQ. 


Figure  4.1:  OTFs  for  three  different  focus  aberrations. 


The  OTF  equation  [15]  for  the  focus  aberration  in  one  dimension  is 


H(fx)  =  A(|^)sinc 


8  W„ 


(4.3) 
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The  effect  of  the  sine  function  on  the  triangle  function  A,  as  the  path  length  error 
Wm  is  increased,  is  clearly  evident  in  the  corresponding  OTF  plots.  fa  is  the  cutoff 
frequency  and  fx  is  the  location  in  the  frequency  domain.  A  is  the  wavelength  of  the 
light.  The  next  section  provides  more  detail  on  the  use  of  the  OTF  in  detection  of 
the  phase  aberration  of  the  EO  image. 

4-2  Focus  Aberration  Detection  Algorithm 

Using  only  post-processing  of  a  camera  image,  the  FAD  algorithm  developed 
estimates  the  amount  of  focus  error  in  the  EO  system.  This  is  accomplished  through 
the  use  of  a  MAP  estimator  [45].  The  FAD  algorithm  detailed  in  this  section  is  shown 
in  Figure  4.2. 

4-2.1  Generate  phase  screens.  The  size  of  the  input  image  is  used  to  deter¬ 
mine  the  size  of  the  aperture  and  the  phase  screen  used.  With  the  assumption  that 
an  image  is  N  by  N  pixels,  the  aperture  and  phase  screen  are  N/2  pixels  in  diame¬ 
ter.  A  Zernike  polynomial  representation  of  the  focus  aberration  phase  screen  is  used 
in  the  ML  estimation  of  the  deconvolved  object.  The  Zernike  polynomial  needs  to 
be  generated  once  at  the  start,  and  then  scaled  by  the  magnitude  of  the  aberration 
throughout  the  remainder  of  the  iterations. 

4-2.2  Input  image.  The  FAD  algorithm  is  initialized  by  reading  in  the 
camera  image  ( d(x,y ))  for  estimation  of  its  focus  aberration.  This  image  can  be  from 
any  arbitrary  scene  under  observation.  An  additional  set  of  dark  images  are  used  to 
set  the  bias  for  the  camera  images  used  in  the  algorithm. 

4-2.3  Generate  log-likelihood  plot.  The  heart  of  the  FAD  algorithm  is  the  use 
of  an  RL  style  algorithm  [30].  The  algorithm  assumes  a  known  aberration.  Starting 
with  a  known  OTF,  the  algorithm  deconvolves  an  object.  The  object  is  then  an  input 
to  a  log-likelihood  equation  to  calculate  a  value  for  that  object.  The  estimator  uses 
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Figure  4.2:  Focus  Aberration  Detection  Algorithm. 


the  PMF  of  a  Poisson  dominated  image 


x  y 


i(x,y)d(x,y^exp  l^x,y^ 
d(x,y)\ 


(4.4) 


where  d  is  a  single  observation  of  a  random  variable  with  range  D  and  z  is  a  single 
observation  of  a  random  variable  with  a  range  Z .  For  this  domain,  d  is  a  single  image 
frame  with  a  range  of  non-negative  integers  and  z  is  the  focus  error  over  the  range  of 
possible  focus  errors.  The  Poisson  PMF  is  used  to  generate  the  log-likelihood  of  the 
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estimate,  C(o,z). 


£(o,z)  =  EE  [d(x,  y)ln(i(x,  y))  -  i(x,  y )  -  ln(d(x,  y)!)]  (4.5) 

x  y 

where  i(x,  y)  is  the  image  generated  from  the  RL  estimated  object  by  convolution  with 
the  known  aberration.  The  term  ln(d(x,y)\)  is  a  constant  term  which  does  not  vary 
as  the  focus  aberration  changes  and  is  therefore  discarded  in  the  calculation.  Due 
to  a  characteristic  of  the  focus  aberration  a  fairly  smooth  plot  of  the  likelihood  of 
the  deconvolved  object  is  generated  as  the  focus  aberration  is  varied.  As  the  known 
aberration  changes,  the  log-likelihood  of  the  estimates  remain  fairly  flat  until  the 
aberration  is  increased  past  the  actual  focus  aberration  present  in  the  system.  Past 
this  point  the  log-likelihood  decreases.  The  algorithm  iterates  in  300  increments  from 
zero  to  some  maximum  focus  aberration.  The  granularity  of  the  setup  used  amounts 
to  208243  waves  or  23nnn  in  displacement  for  each  increment.  This  is  calculated 
using  the  right  side  of  Equation  4.2  with  a  equal  to  the  increment  multiplied  by  the 
atmospheric  Zernike  covariance  value  for  defocus  of  0.0232.  [33].  At  increment  300, 
this  amounts  to  62472814  waves  of  focus  aberration.  The  plots  for  lens  positions  at 
0.439m,  0.388m,  and  0.350m  from  the  CCD  are  shown  in  Figure  4.3(a),  4.3(b),  and 
4.3(c)  respectively. 


4-2-4  Estimate  focus  aberration.  Once  a  log-likelihood  plot  is  generated,  the 
probability  of  the  defocus  estimate  is  added  to  it  using  the  MAP  equation 


Pz\d(Z\D) 


p^(D\Z)P,(Z) 

P*(D) 


(4.6) 


Pz\d(Z\D)  is  the  probability  of  the  focus  aberration  z  given  the  data  d  and  Pz(Z )  is 
a  standard  Gaussian  random  variable  for  the  focus  aberration 


1  -(Z-Z)2 

Pz{Z)  =  — — exp  2^ 
Zna 


(4.7) 
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ML  Poisson 


ML  Poisson 


(a)  0.439m  from  CCD.  (b)  0.388m  from  CCD. 

x-|08  ML  Poisson 


(c)  0.350m  from  CCD. 

Figure  4.3:  Log-likelihood  plots  at  various  lens  positions. 

Z  is  a  random  variable  that  represents  the  actual  focus  aberration  with  a  standard 
deviation  of  a. 

The  log  of  the  MAP  equation  yields 

ln(pAd(Z\D))  =  In  (P^(p|^(Z)) 

=  ln{pdlz(D\Z))+ln(p,(Z))-ln(Pd(D)) 

«  C(o,z)+(Z~J}\  (4.8) 

The  estimate  is  determined  by  an  iterative  algorithm,  depicted  in  the  large  block 
in  the  lower  left  area  of  Figure  4.2,  since  the  actual  value  of  the  mean  focus  aberration, 
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Z,  is  unknown  and  must  be  estimated.  The  maximum  possible  Z  is  used  to  initialize 
the  algorithm.  The  resulting  value  of  Z  that  produces  the  maximum  value  of  the 
sum  is  the  current  estimate  of  the  magnitude  of  the  focus  aberration.  When  the 
current  estimate  differs  from  Z  by  less  than  the  granularity  of  the  plot,  the  algorithm 
terminates.  If  the  difference  is  greater  than  one  increment,  the  estimate  is  then  used 
to  update  the  a  priori  estimate  of  Z  for  the  next  iteration.  For  the  granularity  of 
the  plots  used,  the  algorithm  always  converged  to  a  value  in  less  than  five  iterations. 
This  final  value  is  the  MAP  estimate  of  the  magnitude  of  the  focus  aberration  for  the 
image. 

The  choice  of  value  for  the  focus  covariance  affects  the  performance  of  the  esti¬ 
mator.  If  the  choice  is  too  small  then  poor  estimates  are  possible  as  the  estimator  gets 
trapped  at  a  local  maximum  and  not  the  global  maximum.  The  erroneous  estimates 
are  predominantly  at  the  shorter  focus  aberrations  and  are  seen  in  the  additional 
structure  visible  in  Figure  4.4.  This  is  due  to  the  additional  and  more  pronounced 
local  maximums  in  the  log-likelihood  plots  as  the  actual  focus  aberration  is  decreased. 
An  additional  local  maximum  is  shown  in  Figure  4.3(c)  at  an  aberration  estimate  near 
increment  200.  As  the  choice  of  focus  covariance  is  increased  the  algorithm  will  pass 
over  these  local  maxima  to  find  the  global  maximum.  The  erroneous  estimates  are 
thus  reduced  as  the  focus  covariance  increases  as  shown  in  Figure  4.5.  However,  if 
the  choice  is  too  large,  the  estimator  will  pass  over  all  of  the  values  in  the  likelihood 
plot  and  always  estimate  zero  focus  error.  For  the  remainder  of  the  results,  a  value  of 
212p,m2  is  used  for  the  focus  covariance  as  this  corresponds  to  a  standard  deviation 
of  the  focus  error  of  1.4cm.  This  1.4cm  focus  error  represents  the  uncertainty  in  ones 
ability  to  focus  the  system  manually. 

4-2.5  Repeat  for  remaining  images.  After  the  magnitude  of  the  focus  error 
is  determined,  the  sign  must  be  determined  as  well  in  order  to  move  the  lens  in  the 
correct  direction  to  minimize  it.  In  order  to  accomplish  this  a  second  image  could 
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Lens  distance  to  CCD  (m)  °'35  20  Pixels  across  aperture 

Figure  4.4:  Focus  Aberration  Estimates  versus  Aperture  Size,  Focus  Cov  =  26/rm2. 

be  gathered  after  the  focal  plane  to  lens  relationship  is  altered.  The  resulting  second 
focus  error  estimate  would  reveal  the  sign  of  the  focus  error. 


4-3  Testing  with  Laboratory  Data 

The  laboratory  setup  used  in  the  development  of  the  FAD  algorithm  uses  an 
optics  table  with  a  CCD  camera,  a  compound  lens,  and  a  target  (see  Figure  4.6.)  The 
individual  components  are  described  in  more  detail. 

4-3.1  Target.  A  set  of  bar  targets  provides  an  extended  object  for  testing 
the  algorithm.  A  red  LED  of  700nm  wavelength  is  used  to  backlight  the  bars.  An 
image  of  the  target  is  shown  in  Figure  4.6(a). 

4-3.2  Camera.  The  camera  used  is  a  Cascade  512F  with  a  512  x  512  CCD 
array.  It  is  a  monochromatic  camera  with  a  peak  quantum  efficiency  approximately 
88%  at  700nm.  The  pixel  size  is  16 ym  x  16 ym.  For  the  images  used  in  development 
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Lens  distance  to  CCD  (m) 


0.35 


Pixels  across  aperture 


Figure  4.5:  Focus  Aberration  Estimates  versus  Aperture  Size,  Focus  Cov  =  212 /im2. 


of  the  defocus  detection  algorithm,  the  integration  time  is  set  at  1/zs  and  only  the 
central  256X256  pixels  are  used.  The  resulting  images  have  a  signal-to-noise  (SNR) 
ratio  for  the  various  positions  shown  in  Figure  4.7. 


4-3.3  Lens.  The  lens,  shown  in  Figure  4.6(b),  is  composed  of  two  convex 
lenses  with  a  focal  length,  /,  of  400 mm.  The  lenses  are  mounted  on  a  translation 
table  as  close  together  as  possible.  The  separation  between  the  two  lenses,  dierLS,  is 
approximately  31mm.  An  aperture  of  4mm  in  diameter  is  placed  in  front  of  the  lens 
assembly.  The  compound  focal  length  is  approximately  192mm  using  the  equation 
for  the  front  focal  length  [16] 


f-f-l. 


f Miens  ~  /2) 

dlens  (/i  +  A) 


(4.9) 


of  a  compound  lens.  Because  the  focal  length  of  both  lenses  is  the  same,  f\  —  j'2  —  /, 
both  the  front  focal  length  and  the  back  focal  length  are  the  same. 
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(a)  Target.  (b)  Compound  Lens.  (c)  Camera. 


Figure  4.6:  Laboratory  Components. 

4.3.4  Images.  A  series  of  images  are  collected  with  the  compound  lens 
adjusted  in  13 mm  increments  from  an  arbitrary  starting  location  at  0.45m  in  front 
of  the  camera.  A  set  of  images  is  taken  at  that  location  and  then  the  lens  is  moved 
towards  the  camera.  Another  set  of  images  is  then  taken  at  the  new  location.  This  is 
repeated  for  an  arbitrary  number  of  locations.  Only  the  compound  lens  is  mounted  to 
a  translation  table.  Both  the  target  and  the  camera  are  mounted  in  fixed  locations. 


4-4  Additional  Considerations 

4-4-1  Sampling  effects.  As  with  any  image  processing  algorithm  sampling 
has  an  effect  on  the  results.  To  determine  the  minimum  sampling  requirements  in  the 
digital  domain,  both  the  cutoff  frequency  of  the  camera  CCD  array  and  the  aperture 
are  analyzed.  First,  the  cutoff  frequency  of  the  camera  is  fixed  due  to  the  pixel  pitch 
(16  fimX  16 /im)  of  the  CCD  array.  This  fixes  the  sampling  frequency  of  the  camera  at 

CCDfs  =  •  10e6  m”1  =  62500m”1  (4.10) 


for  the  CCD.  The  images  used  are  256  x  256  pixels.  Dividing  by  the  number  of 
pixels  in  the  image  (. N  =  256)  results  in  the  contribution  to  the  spectrum  per  pixel, 
Sp  =  244m”1.  Second,  the  cutoff  frequency  for  incoherent  imaging  is  defined  as  [15] 
fc  =  2  f0  where 


fo 


D 

2X-z' 


(4.11) 
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Figure  4.7:  Signal-to- noise  ratio  versus  lens  position. 

D  is  the  aperture  diameter  and  z  is  the  distance  between  the  CCD  and  the  aperture. 
With  the  aperture  at  a  distance  of  0.45m  from  the  CCD,  the  cutoff  frequency  is 
12855m-1.  With  the  aperture  at  a  distance  of  0.33m  from  the  CCD,  the  cutoff 
frequency  is  17305m-1.  Thus  the  constraint  on  the  cutoff  frequency  of  the  image  is 
determined  by  the  location  of  the  lens  array. 

To  meet  minimum  sampling  requirements  in  the  digital  domain,  the  number  of 
samples  across  the  “aperture”  is  calculated  by  dividing  the  limiting  cutoff  frequency 
by  the  spectrum  per  pixel  contribution,  Sp.  For  the  aperture  at  its  farthest  point  in 
the  laboratory  setup,  a  minimum  of  52  pixels  across  the  digital  aperture  is  required. 
At  its  nearest  point,  a  minimum  of  70  pixels  is  required.  The  effects  of  sampling  on 
the  corresponding  estimates  are  shown  in  Figure  4.5.  The  effects  of  undersampling 
are  most  noticeable  for  positions  1  to  4  if  the  aperture  is  less  than  50  pixels  across. 

To  explore  the  effects  of  sampling  the  number  of  pixels  across  the  digital  aperture 
is  varied  from  16  to  128  pixels.  An  interesting  effect  occurs  with  the  focus  aberrations 
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(a)  Position  1  (,45m) 


(d)  Position  4  (,41m) 


(g)  Position  7  (.38m) 


(b)  Position  2  (.44m) 


(e)  Position  5  (,40m) 


000 


(h)  Position  8  (.36m) 


(c)  Position  3  (.43m) 


(f)  Position  6  (.39m) 


D0O 


(i)  Position  9  (.35m) 


Figure  4.8:  Camera  images  at  various  positions  (distance  to  CCD.) 


estimates  as  the  number  of  pixels  across  the  aperture  is  changed.  Although  the  focus 
aberration  estimates  look  like  they  are  increasing  as  the  number  of  aperture  pixels 
are  increased,  the  net  effects  on  the  deconvolved  objects  are  negligible.  Looking  at 
Equation  4.3  again,  this  is  explained  by  the  change  in  fx  as  the  number  of  pixels 
is  changed.  As  the  number  of  pixels  in  the  digital  “aperture”  is  increased  there  is  a 
corresponding  decrease  in  frequency  steps,  fx,  in  the  frequency  domain.  Since  fa  is 
fixed  and  fx  changes  with  the  number  of  pixels,  a  corresponding  change  in  Wm  is 
required  to  maintain  the  same  OTF  as  the  sampling  changes. 

This  is  demonstrated  by  looking  at  the  difference  in  deconvolved  objects  with 
an  aperture  width  of  70  pixels  versus  an  aperture  width  of  128  pixels.  First,  setting 
the  width  of  the  aperture  to  70,  as  predicted  by  the  previous  analysis,  results  in  the 
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(a)  Position  1  (b)  Position  3  (c)  Position  5 


(d)  Position  7  (e)  Position  9 


Figure  4.9:  Objects  using  70  pixel  wide  aperture. 

corresponding  deconvolutions  shown  in  Figure  4.9.  Changing  the  aperture  width  to 
128  pixels  across  results  in  the  deconvolutions  shown  in  Figure  4.10.  One  can  see  that 
there  is  little  difference  between  the  sets  of  deconvolved  objects.  The  only  difference 
visible  to  the  eye  is  between  the  two  objects  at  position  5. 


(a)  Position  1  (b)  Position  3  (c)  Position  5 


(d)  Position  7  (e)  Position  9 

Figure  4.10:  Objects  using  128  pixel  wide  aperture. 


If  one  takes  estimates  that  were  made  at  one  aperture  width  and  then  uses 
them  to  deconvolve  objects  using  an  aperture  with  a  different  width,  the  results 
are  visually  poor.  The  larger  the  actual  aberration  and  the  greater  the  difference 
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in  aperture  width,  the  greater  the  error  in  deconvolution.  This  is  demonstrated  by 
taking  the  estimates  from  the  FAD  algorithm  with  an  aperture  width  of  70  pixels  and 
using  them  to  deconvolve  objects  using  an  aperture  width  of  128  pixels.  The  resulting 
objects  are  shown  in  Figure  4.11.  The  most  significant  errors  are  at  positions  1  and 
3  which  are  at  the  greatest  distance  from  the  lens  and  thus  have  the  largest  actual 
focus  aberration  in  the  images. 


t  i  I  ® 


(a)  Position  1  (b)  Position  3  (c)  Position  5 


010  I  III 


(d)  Position  7  (e)  Position  9 

Figure  4.11:  Objects  using  70  pixel  estimates  with  128  pixel  wide  aperture. 

The  results  of  looking  at  the  sampling  requirements  are  what  one  expects  from 
sampling  theory.  The  first  result  is  that  undersampling  results  in  poor  estimates. 
The  second  result  is  that  once  sufficient  sampling,  i.e.  the  Nyquist  frequency  [15],  is 
reached  additional  oversampling  has  little  effect  on  the  results. 

4-4-2  Autofocus  Algorithm.  With  at  least  one  more  image  taken  at  a  differ¬ 
ent  lens  location,  the  focus  aberration  estimates  can  be  used  to  determine  which  side 
of  the  focal  point  the  lens  is  on  and  thus  the  direction  and  relative  movement  of  the 
lens  required  to  focus  the  image.  Using  the  estimated  aberration  for  each  image,  the 
lens  can  then  be  repositioned  to  minimize  the  focus  aberration. 

4-4-3  Speedup.  Several  considerations  have  an  effect  on  the  wall-clock  time 
required  by  the  FAD  algorithm.  As  written  for  this  research,  for  each  image  processed, 


83 


the  FAD  algorithm  requires  approximately  75  minutes  of  wall-clock  time  running  on 
a  64-bit  processor  running  at  3GHz.  The  breakdown  of  the  timing  primarily  comes 
from  loading  the  images  which  requires  2.5  seconds,  15  seconds  per  iteration  of  the 
log-likelihood  plot  generation  loop,  less  than  a  second  for  the  rest  of  the  algorithm. 
Therefore,  the  most  significant  way  to  speed  up  the  algorithm  is  to  speed  up  the  log- 
likelihood  plot  generation  loop.  The  number  of  loops  required  by  the  RL  algorithm  is 
at  a  minimum.  Any  additional  reduction  of  the  stall  condition  for  the  RL  algorithm 
introduces  additional  structure  in  the  plots  that  cause  the  FAD  algorithm  to  produce 
erroneous  estimations.  A  faster  alternative  would  be  to  explore  the  use  of  an  inverse 
Weiner  filter  to  deconvolve  the  object  instead  of  the  RL  algorithm. 

There  are  several  other  possibilities  to  decrease  the  algorithm’s  wall-clock  time 
depending  on  application.  One  simple  improvement  for  immediate  speedup  comes 
from  a  priori  knowledge  of  the  focus  aberration  limit  inherent  in  the  electro-optical 
system.  For  the  current  setup  this  would  result  in  reducing  the  plot  generation  loop 
from  300  to  200  points.  That  would  immediately  reduce  the  wall-clock  time  by  one- 
third  to  approximately  50  minutes  for  the  laboratory  setup  used.  Another  improve¬ 
ment  by  doubling  the  granularity  of  the  plots  would  halve  the  wall-clock  time.  The 
trade  space  in  the  algorithm  can  be  tailored  to  the  application.  One  additional  mod¬ 
ification  is  valid  regardless  of  application.  Near  linear  speed  up  can  be  achieved  with 
additional  processors  used  to  produce  the  ML  plot.  Adding  one  additional  processor 
cuts  the  time  required  to  produce  the  plot  basically  in  half.  The  limit  is  attained 
when  a  processor  is  available  for  every  point  in  the  ML  plot.  This  would  result  in 
a  wall-clock  time  of  approximately  18  seconds  under  the  current  setup.  Additional 
speedup  can  be  attained  by  a  faster  processor  or  specialized  hardware  tailored  for  the 
algorithm. 

This  chapter  discussed  the  production  of  test  data  with  calibrated  focus  errors. 
To  detect  the  focus  aberration  present  in  the  laboratory  images,  the  FAD  algorithm 
was  developed.  The  calibrated  data  is  available  for  testing  of  the  algorithms  derived 
in  the  next  chapter. 
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V.  Polarimeter  Image  Diversity  Blind  Deconvolution 
Algorithm  Development 


Using  historical  development  of  blind  deconvolution,  presented  in  Section  2.2, 
a  blind  deconvolution  algorithm  for  a  two  channel  polarimeter  is  developed. 
The  algorithm  makes  use  of  an  additional  constraint  between  two  imaging  channels. 
The  general  development  follows  the  generalized  expectation  maximization  (GEM) 
algorithm  development  by  Schulz  [36]  for  his  multiframe  blind  deconvolution  (MFBD). 

Starting  with  the  historical  basis  for  a  single  channel,  the  statistical  model  is 
extended  to  a  two-channel  polarimeter  under  the  assumption  of  a  known  psf.  This 
is  followed  by  development  of  a  polarimeter  MFBD  (PMFBD)  algorithm  when  the 
assumption  is  removed. 


5.1  Historical  Basis 


Standard  imaging  systems  have  a  PMF  for  a  single  image  such  as 


pmf(c) -nil 


x  y 


,i(x,  y)d(x,y^exp  l(x’y\ 
d(x,  y)\ 


(5.1) 


for  Poisson  distributed  intensity  images.  Where 


i{x,  =  °(Z1  w)h(x  -  z,y-w).  (5.2) 

2  W 

Some  terms  defined  for  use  in  the  GEM  algorithm  developed  are: 

•  {d(x,y\z,w)}  is  the  complete  data  (CD)  set  and 

•  {d(x,  y)}  is  the  incomplete  data  set 

where  there  is  a  many-to-one  mapping  between  the  two  sets  as  defined  by  d(x,y)  = 
x,y\z,w).  Additionally,  the  expectations  of  the  data  are  defined  by  Schulz 
as 


E[d(x,y\z,w)\  =  o(z,w)h(x  —  z,y  —  w) 


(5.3) 


85 


and 


E[d(x,y)}  =  ^2^2E[d(x,y\z,w)] 

Z  W 

=  ^2^2°(z,w)h(x  -  z,y  -w) 


Z  W 

=  i(x,y,o). 

The  CD  log-likelihood  function  of  the  PMF  is 

L CD{o)  =  ^2^2^2^2[d(x,y\z,w)ln(o{z,w)h(x  -  z,y -w)) 

z  w  x  y 

—o(z,  w)h(x  —  z,y  —  w)  —  ln(d(a;,  y\z,  to)!)]. 


(5.4) 


(5.5) 


The  last  term  is  a  constant  and  is  therefore  dropped  when  maximizing  the  log- 
likelihood  equation  becoming 

L CD(o)  =  | d(x,  y\z ,  w)ln(o(z ,  w)h(x  —  z,y  —  w))  —  0(2,  w)h{x  —  z,y  —  w)  .  (5.6) 

z  w  x  y 

A  GEM  algorithm  for  deconvolution  with  a  known  psf,  requires  that 


LCD(()ne  w)  >  jCD  ,old 5 


(5.7) 


and 


where 


Q(onew\o°ld )  >  Q(o°ld\o°ld) 


old  I  ~old\ 


Q(o\o°ld)  =  Eold  [LCD(o)\{d(x,y)}' 


Maximizing  Q  with  respect  to  the  object  results  in 


(5.8) 


(5.9) 


</"”  -  urn  in;:)'.: (!\o  o°“). 
oeD 


(5.10) 


To  satisfy  the  GEM  requirement,  one  takes  the  partial  derivative  of  the  complete 
log-likelihood  function  with  respect  to  the  object  and  maximizes  it.  The  partial 


derivative  is 


dLCD 
d o(z0,  w0) 
d 


do(z0,w0 


EEEE  [d(x,  y\z,  w)ln(o(z ,  w)h(x  —  z,y  —  w)) 


z  w  x  y 

—o(z,  w)h(x  —  z,y  —  w)] 


EE 

X  y  L 

EE 

x  y 

EE 


rf(x,  r/|z,  w)  d 
o(z0,w0)  do(z0,w0 ) 

_ 

(9o(^0,wc 

_ 

L  o(z0,w0) 

d(x,y\z,w ) 
o(^0,w0) 


x  1/ 


X!  o(z,  w)h(x  -  z,  y  - 


h(x  -za,y-  w0) 


Y1  ~z^y~  w°)\ 


w 


x  y 


o(za,  W0) 


-  1 


(5.11) 


The  result  of  the  partial  derivative  is  set  equal  to  zero  and  solved  for  o: 


o(z0,  w0)  =  ^  ^  d(x,  y\z ,  w).  (5.12) 

x  y 

Taking  the  expectation  of  the  result  is  the  expectation  step  in  the  EM  algorithm.  It 
is  then  maximized  using  the  results  of  Shepp  and  Vardi.  The  expectation  given  the 
old  object  is 

noidiM  I  MM  M  h(x-z0,y-w0)oold(x,y) 

E  [d(x,y\z,w)\d(x,y)]  = - - - -gz - d(x,y).  (5.13) 

t(x,y,oOLa) 

Substituting  this  into  Equation  5.12  yields 


o 


new 


d(x,y) 
i(x,y;  oold) 


h(x  —  z0,  y  —  w0). 


which  is  the  RL  algorithm. 


(5.14) 
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5.2  Algorithm  Development 


The  first  step  in  developing  a  polarimeter  blind  deconvolution  algorithm  is  to 
develop  a  polarimeter  deconvolution  algorithm  using  a  known  psf.  The  approach 
taken  is  to  modify  the  conventional  GEM  algorithm  to  work  with  two-channel  po¬ 
larimeter  data.  The  PMF  for  the  combination  of  a  conventional  image  with  an  added 
polarization  image  is 


PMF{0.p) 


nnnn 

z  w  x  y 


fih(x,y;  o,p,  h)dh(x,y^exp  lh{x.y\o,P,h) 
dh(x,y)\ 

i(x ,  y\  o ,  hy(x’y) exp~dx’y'°’} d , 
d{x,y)\ 


(5.15) 


The  image  definitions  for  notation  purposes  arbitrarily  define  one  image  as  the  hor¬ 
izontally  (h)  polarized  image  and  the  other  as  the  unpolarized  primary  image.  The 
individual  images  are  defined  as 


i(x,y,o)  =  ih(x,y,o,p)  +iv(x,y,o,p) 

=  7}  ^2  °(z°’  w° )  (1  +  Wo))hi  (x  -  z0,y-  wa) 

x  y 

+  ^  ^2^2o(z0,w0)(l  -  V(z0,w0))h2(x  -  za,y  -  wa))  (5.16) 

x  y 

where  h(x,y )  is  the  psf  and  V(z,w)  is  the  degree  of  polarization  [14]. 

5.2.1  Single  polarimeter  channel  with  known  psf.  The  log-likelihood  equa¬ 
tion  for  a  single  polarimeter  channel  combined  with  the  conventional  channel  is 

LOD(o,p) 

=  ££££  [dh(x,  y\z,  w)ln(ih(x,  y))  -  ih(x,  y) 

z  w  x  y 

+d{x,  y\z,  w)ln{i{x,  y))  -  i{x,  y)\ 
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(5.17) 


=  E  E E  E^1’  y\z’  w)ln(o(z,  w )p(z,  m)h(x  -  z,  y  -  m)) 

z  w  x  y 

—o(z,  w)p(z,  w)h(x  —  z,y  —  w) 

+d(x,  y\z,  w)ln(o(z,  w)h(x  —  z,y  —  w)) 

—o(z,  w)h(x  —  z,y  —  w)] 


where  p(z,w)  =  1  +V(z,  w)  to  simplify  the  derivation.  Taking  the  partial  derivative 
with  respect  to  o  yields 


dL  (o,p) 
do(z0,  w0) 
d 


do(z0,w0 


z  w  x  y 


EEEEw*  ,  y\z,  w)ln(o(z ,  w)p(z,  w)h(x  —  z,y  —  w)) 

i 

—o(z,  w)p(z,  w)h(x  —  z,y  —  w) 

+d(x,  y\z,  w)ln(o(z,  w)h(x  —  z,y  —  w )) 


EE 


x  y 


EE 


x  y 


dh(x,  y\z,  w)  _  _ 

L  o(z0,w0)  {do(z0,w0 ) 

d 


o(z,  w)h(x  —  z,y  —  w)] 

9  (  \\ 
o[z,w)) 


do(z0,  wa 
d(x,y\z,w ) 


o(z,  w)p(z,  w)h(x  —  z,y  —  w) 
d 


L  o(^o,  Wo)  <9o(^0,  Wo) 

d 


o(z,w)) 


EE 


a;  y 


do(z0,w0) 

dh(x,y\z,w 


[  o(z0,w0) 


o(z ,  w)h(x  —  z,y  —  w) 

-  ~  p{z0,  w0)h(x  -  zQ,y~  w0) 


d(x,y\z,w )  ' 

H - 7 - 7 - h(x  -Zo,y-  Wo) 


o(z0,  Wo) 

J2x  4(g,  2/|z,  w)  +  E.r  Ey  d(x,  y\z,  w) 

o(z0,  Wo) 


-p(z0,w0)  -  1. 


(5.18) 


Setting  this  equal  to  zero  and  solving  for  o(z,  w ) 


0  = 


Ex  Ej,  40,  y\z,  w )  +  Ex  E„  d(x,  y\z ,  w) 


x  ^ y 


1  +  p(z0,  W0) 
0(za,  W0) 


o(z0,  Wo) 

Ex  Ey  40,  j/Q  4  +  Ex  Ey  d(x,  y\z,  w) 

o(z0,  Wo) 

Ex  Ey  40,  4b  4  +  Ex  Ey  <*P,  s/I*,  4 

1  +  p(z0,w0) 


p(z0,w0 )  -  1 


(5.19) 


Taking  the  expectation  given  oo  d  and  pod  yields 


nold  .  V  V  dh{x,y)pold(z,w)u/  „ 

0  Z^X  2^2/  ih(x,y,oold,pold)  'H3'  Z,y  w / 

1  +  pp,  iu) 


=  o 


o/d 


1  +  pp,  it/ 

^  f  pold(z,w)-dh(x,y) 

L-^dX  E—/1 


y  l  ih(x,y;oold,pold)  i(x,y,oold) 


+ 


d(x,y) 


h(x  —  z,y  —  w) 


1  +  pp,  in) 


(5.20) 


In  order  to  implement  this  an  estimate  for  p(z,  w)  in  the  denominator  must  be 
determined.  This  causes  a  slight  deviation  from  a  simple  GEM  algorithm.  In  this 
case,  p(z,w)  is  set  to  pold  with  good  results  as  shown  in  Section  6.1.2.  The  iterative 
algorithm  becomes 


^ new  _  qOIcL 


Ex  E, 


ih  (x,y,oold,pold) 


+ 


d(x,y) 

i(x,y;oold) 


h(x 


z,y 


w 


1  +  pold( S 


(5.21) 
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Similarly,  taking  the  partial  derivative  with  respect  to  p  yields 


<9L  (o,p) 
dp(z0,w0 ) 
d 


dp(zQ ,  wQ 


EE 

a;  y 

d 


dp(zQ,w , 


■EEw*  ,  y|z,  w)ln(o(z ,  w)p(z ,  w)h(x  —  z,y  —  w )) 

x  3/ 

— 0(2;,  w)p(z,  w)h(x  —  z,y  —  w) 

+d(x,  y\z,  w)ln(o(z ,  w)h(x  —  z,y  —  w)) 
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EE 
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x  y  L 

Ex 


o(z0,  w0)h(x  -  zQ,y  -  wa) 


p(z0 ,  w0 


-  o(20,iy0). 


(5.22) 


Setting  this  equal  to  zero  and  solving  for  p(z,  w ) 


p{z0,w0)  = 


Ex  Eyjfe^ky^) 
0(^0,  Wo) 


(5.23) 


Using  the  expectation  of  dh(x,y\z,w)  given  oold  and  pold  yields 


P ' 


Ex  Ey  [4(g,  y|u  w)  |4(u  y)\ 

o(z0,w0) 

h(x-Zo,y-Wo)0old(Zo,Wo)pold(Zo,Wo)  1  /  \ 

it  (x ,y,o°^d ,p°^)  5  U J 


poWEE 

x  3/ 


0(^0,  WD) 

0old(zo,wo)dh(x,y ) 
o(z0,w0)ih(x,y ;  0old,pold' 


■h(x  —  z0,y  —  w0)  (5.24) 
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with  an  assumption  similar  to  that  used  in  Equation  5.21,  such  that  o(z0,w0 )  in  the 
denominator  is  set  to  oold,  this  simplifies  to 


piiew  _  p°ld 


x  y 


dh(x,y) 


ih(x,  y,  oold,pold 


■h(x  -  Zo,y-  w0). 


(5.25) 


Using  the  results  of  Equations  5.18  and  5.25  a  new  GEM  algorithm  is  created  for 
polarimeter  deconvolution  using  a  known  psf. 


5.2.2  PMFBD  GEM  algorithm.  The  polarimeter  deconvolution  algorithm  is 
extended  to  blind  deconvolution  using  a  methodology  similar  to  Schulz  by  adding  mul¬ 
tiple  frames.  The  log-likelihood  equations  developed  in  the  previous  section  are  used 
to  develop  a  multiframe  GEM  algorithm  for  blind  deconvolution  using  polarimeter 
information.  The  desired  property  is  that 


L (0new,pnew,hnew)  >  L (oold,pold,hold) 


(5.26) 


and 


Q(o,p,  h\o°u,p°,d,  hM)  =  Qo(o\o0,d,p0,\  h°ld)+Qp(j>\o°ld,p°ld,  h°,d)+Qh(h\o°'d,p°'d,  h°,d) 

(5.27) 

where 


Q(o,p,h\om,pM,hm)  =  Em  [lPD(o,p,h)\{drt(x.y)}.{Ak(x,y)}Y  (5.28) 

where  Eold  [-\{dpk(x,  y ),  dk(x,  y)}]  is  the  expectation,  conditioned  on  the  data  from  one 
polarization  channel  dpk(x,y )  and  the  data  from  the  primary  channel  dk(x,y),  and 
assuming  that  the  underlying  object  is  oold  and  the  polarization  data  is  pold.  Each 
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image  pair,  dpn  and  dm  are  assumed  to  have  the  same  psf  if  n  =  m  and  different  psfs 
if  n  ^  m. 

Any  algorithm  that  meets  these  requirements  such  that 


Q0(onew\oold,p°ld,h°ld) 
Qp(pnew\oold,pold,hold) 
Qh(hnew\oold,pold,  hold) 

where 


>  Q0(oold\oold,pold,h°ld )  (5.29) 

>  Qp(pold\oold,pold,hold)  (5.30) 

>  Qh(hold\oold,p°ld,  hold)  (5.31) 


0 new 

=  arg  max  Q0(o\o°ld ,  pold ,  hold ) 

oeO 

(5.32) 

pn€W 

=  arg  max  Qp(p\o°ld,  p°ld,  hold) 

(5.33) 

ypnew 

=  arg  max  Qh  (h\oold,  pold,  hold ). 
hesj 

(5.34) 

is  a  GEM  algorithm. 

The  complete  log-likelihood  algorithm  for  one  polarization  channel  in  combina¬ 
tion  with  the  primary  channel  with  multiple  frames  is 


L rD{o,p,  hk) 
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(5.35) 
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The  previous  results  for  the  deconvolution  with  a  known  psf  still  hold  for  cal¬ 
culating  onew  and  pneu\  there  are  just  individual  loops  for  each  frame  and  a  method 
for  combining  the  results.  In  this  case  a  simple  averaging  is  used  such  that 


mean(o^ew) 


pnew  =  mean(pnewj 


(5.36) 

(5.37) 


However,  an  iterative  update  is  derived  for  the  individual  psfs  by  maximizing  the 
complete  log-likelihood  equation  with  respect  to  the  psfs.  The  partial  derivative  of 
the  complete  log- likelihood  equation  with  respect  to  hk(z0,w0 )  is 


d 


dhk(x0,y0) 

d 

dhk(x0,y0) 


L  (o,p,  hk)  = 

£££££  dpkfai  y\-£  %->y 

L  k  z  w  x  y 

■  ln(o(x  —  z,y  —  w)p(x  —  z,y  —  w)hk{z ,  w)) 

°(x  ~ z' y  ~  w^x  ~  z> y  ~  w) 

k  z  w  x  y 

£££££  dk(x,  y\x  -  z,y  -  w)  ln(o(x  -  z,y  ~  w)hk(z,  w)) 


k  z  w  x  y 


££ 


-  °(x  ~ z’ y  ~  w)h^  w) 
k  z  w  x  y 

dpk(x,y\x  -  z,y-w) 


hk  (*^05  Vo 


(5.38) 

-  y,  y,  o(x0  —  z,yQ  —  w)p(x0  -  z,yQ~  w ) 

(5.39) 


dk(x,  y\x  —  z,  y  —  w)  /  m 

hk{x0,y0) 


Setting  this  equal  to  zero  yields 


hk  (x0)  yQ ) 


E»  dpk(x,  y\x  -  z,y  -w)  +  Ez  Ew  dk(x,  yfx  -z,y-w) 
E *  E„,  °(xo  ~z,y0-  w) (1  +  p{x0  -  z,yQ-w )) 


(5.40) 
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Taking  the  expectation  given  oold,  pold ,  and  hfd  yields 


Kew&y)  = 


EE 


h1ld(z,w)ooldlx—z,y—w)pold(x—z,y—w)  ,  /  \ 

ipk(x,y;0°ld,poldthold)  dpk(X,  y) 

o{z,w)(  1  +p(z,w)) 

h?ld(z,w)oold(x—z,y—w)  7  /  \ 

7,  (V  «rnold  uold\  d}z\X)  y) 


EE 

2  w 


ik(x,y;oold,hold) 

o(z,w)(  1  +  p(z,  w)) 


(5.41) 


With  the  same  assumptions  on  o(z,w)  and  p(z,w )  as  in  the  previous  development 
this  simplifies  to 


hnew 


tf'-EE 


=  ftf-EE 


dpk(x,y)pold{x—z,y—w)oold(x—z,y—w)  ^  dk(x,y)oold(x—z,y—w) 


ipk  (x  ,y;oold  ,pold  ,hold) 


ik(x,y,oold,hold) 


oold(z,w)(  1  +poW(z,w)) 

-z,y-  w)oold(x  -  z,y-w) 


ipk  (x  ,y;oold  ,pold  ,hold)  * 


oold(z,  ty)(l  +  pold(z,  w )) 


dk (x.y) _  oldf^,  ^  _  ri.A 

ik(x,y;oold,hold)°  \X  V  W) 

oold(z,w)(  1  +  poW(z,  w)) 


(5.42) 


This  chapter  developed  two  polarimeter  deconvolution  algorithms  for  a  two- 
channel  polarimeter,  where  one  channel  is  unpolarized  and  the  other  is  polarization 
sensitive  with  arbitrary  orientation.  The  initial  algorithm  is  an  extension  of  the  RL 
algorithm  to  include  the  additional  polarimeter  channel  assuming  a  known  psf.  The 
second  is  based  on  a  GEM  approach  and  is  a  polarimeter  multiframe  blind  decon¬ 
volution  algorithm.  The  previously  developed  calibrated  data  is  used  to  test  the 
algorithms.  The  results  are  presented  in  the  next  chapter. 
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VI.  Results  and  Analysis  of  Polarimeter  Algorithms 

The  results  and  analysis  of  the  algorithm  development  in  Chapter  V  are  presented 
in  this  chapter.  Comparison  of  the  RL  algorithm  with  the  polarimeter  deconvo¬ 
lution  algorithm  using  simulation  is  presented  in  the  first  section.  This  is  followed  by 
the  results  using  actual  laboratory  images.  The  final  section  presents  a  comparison 
between  the  Polarimeter  MFBD  algorithm  developed  and  an  RL  based  MFBD. 

In  order  to  compare  the  two  algorithms’  performance  a  metric  is  required.  The 
metric  chosen  is  based  on  the  fidelity  metric  of  Kattnig  et  al.  [22],  To  compare  the 
two  algorithms  it  is  implemented  as 


/  ,  •  ,  t  •  ,  \  _  ||  object  obj  CCt  deconvolved,  \  \  L2  Sc  1\ 

J  xddzty yOOJ CCt ,  00 J CCt deconvolved)  Tj  j  "  7779  • 

\\°b]ect\\l2 

A  perfectly  reconstructed  object  would  result  in  a  value  of  zero  for  the  metric. 
Therefore  when  comparing  the  results  of  both  algorithms,  a  lower  value  is  a  better 
deconvolution.  This  applies  only  to  the  simulation  results  since  the  true  object  is 
unknown  when  using  the  laboratory  images.  Using  simulated  images  that  only  have 
diffraction  present,  the  average  of  100  images  results  in  the  diffraction  limited  fidelity 
scores  shown  in  Figure  6.1  for  digital  aperture  sizes  of  32,  48,  and  64  pixels  across 
the  simulated  aperture.  The  images  used  are  128  x  128  pixels  in  size.  The  results 
show  that  as  the  aperture  size  is  decreased  the  resulting  image  is  diffracted  more  thus 
resulting  in  a  lower  fidelity  score. 

6.1  Known  psf  with  simulation  data 

A  Matlab®  simulation  is  used  to  compare  the  new  polarization  sensitive  algo¬ 
rithm  with  the  traditional  RL  algorithm.  The  simulated  target  is  a  bar  target  with 
a  variable  separation  between  the  two  bars.  The  bars  are  of  variable  widths  in  pixels 
and  a  fixed  17  pixels  in  height.  The  fidelity  mapping  algorithm  starts  at  a  minimum 
separation  of  one  pixel  and  increments  up  to  a  16-pixcl  separation.  The  algorithm 
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Figure  6.1:  Fidelity  metric  for  diffraction  limited  images  for  various  digital  aperture 
sizes,  using  target  with  bar  width  of  3  pixels,  target  irradiance  of  1000 
photons  per  pixel,  and  image  size  of  128. 


also  increments  through  200  steps  of  increasing  focus  aberration.  This  is  repeated  for 
various  noise  levels.  The  SNR  [9]  for  the  test  images  is  calculated  using  the  equation 


SNR 


mean(signal) 
standard  deviation(noise ) 

2  *  bar  width  *  bar  height  *  number  of  photons  per  pixel 
standard  deviation(noise)  *  N2 


(6.2) 


The  standard  deviation  of  the  noise  for  the  Matlab®  simulations  is  fixed  at  0.6.  Figure 
6.2  shows  the  results  for  a  bar  width  of  three  pixels,  an  irradiance  of  1000  photons  per 
pixel,  and  an  image  size  of  128  pixels.  The  SNR  of  the  test  images  is  approximately 
10.  In  the  simulations,  the  SNR  is  calculated  from  the  RL  image,  the  other  two  images 
have  half  the  value.  This  simulates  the  use  of  a  50/50  beam  splitter  cube  in  the  optical 
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(a)  RL  fidelity.  (b)  Polarimeter  fidelity.  (c)  RL  fidelity  -  Polarimeter  fi¬ 

delity. 


Figure  6.2:  Map  of  fidelity  metric  for  both  the  RL  algorithm,  the  polarimeter  algo¬ 
rithm,  and  their  difference  at  an  SNR  of  approximately  10. 

path.  With  the  configuration  used  this  results  in  the  sum  of  the  polarimeter  images 
having  3/4  of  the  photons  of  the  RL  algorithm. 


(a)  RL  fidelity.  (b)  Polarimeter  fidelity.  (c)  RL  fidelity  -  Polarimeter  fi¬ 

delity. 


Figure  6.3:  Map  of  fidelity  metric  for  both  the  RL  algorithm,  the  polarimeter  algo¬ 
rithm,  and  their  difference  at  an  SNR  of  approximately  1. 

Changing  the  signal  or  the  irradiance  of  the  target  causes  a  corresponding  change 
in  the  effectiveness  of  both  algorithms.  Reducing  the  signal  to  a  target  irradiance 
of  100  photons  per  pixel  reduces  the  SNR  of  the  images  to  approximately  1.  The 
changes  in  the  fidelity  metric  are  shown  in  Figure  6.3.  As  the  aberration  is  increased 
there  is  a  corresponding  decrease  in  performance  of  both  algorithms.  The  greater 
the  aberration,  the  more  spread  out  the  photons  are  in  the  resulting  image  and  thus 
harder  to  reconstruct  the  original  object  due  to  a  lower  image  SNR.  If  the  assumption 
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of  a  beamsplitter  is  removed  and  one  is  considering  the  difference  between  the  two 
algorithms  on  the  same  image,  the  results  are  shown  in  Figure  6.4.  Note  that  in 
Figure  6.4(a)  the  RL  algorithm  now  has  half  of  the  photons  of  the  previous  results 
shown  in  Figure  6.3(a).  This  results  in  the  failure  of  the  RL  algorithm  at  much  lower 
aberrations  than  the  polarimeter  algorithm.  Also,  instead  of  the  RL  algorithm  having 
better  performance  in  the  low  aberration  region  due  to  more  light,  the  polarimeter 
algorithm  now  performs  better  than  the  RL  algorithm  in  this  region.  It  is  clear  that 
the  polarimeter  algorithm  performs  better  than  the  RL  algorithm  when  working  with 
the  same  image  from  the  primary  channel. 


(a)  RL  fidelity.  (b)  Polarimeter  fidelity.  (c)  RL  fidelity  -  Polarimeter  fi¬ 

delity. 

Figure  6.4:  Map  of  fidelity  metric  for  both  algorithms  without  beamsplitter  assump¬ 
tion  at  an  SNR  of  approximately  1. 

As  the  light  levels  are  reduced,  the  polarimeter  algorithm  continues  to  maintain 
a  better  fidelity  than  the  RL  algorithm.  Reducing  the  number  of  photons  per  pixel 
of  the  target  to  only  10  photons  per  pixel  results  in  an  SNR  of  approximately  0.1. 
The  cross  sectional  views  of  the  deconvolved  objects  by  both  algorithms,  in  Figure 
6.5  shows  an  interesting  difference  between  the  two  algorithms.  For  this  comparison 
polarimeter  algorithm  has  the  beamsplitter  assumption  applied  and  even  with  only 
3/4  of  the  photons  that  the  RL  algorithm  has,  the  cross  sections  clearly  show  a  much 
greater  contrast  for  the  polarimeter  algorithm. 
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Figure  6.5:  Cross-sectional  views  of  deconvolved  object  with  a  barwidth  of  3  pixels 
and  a  target  irradiance  of  10  photons  per  pixel. 

6.1.1  Comparison  of  results  with  CRLB  predictions.  The  bounds  calculated 
in  Section  3.3.11  show  the  ability  of  the  polarimeter  algorithm  to  produce  a  better 
estimate  than  conventional  imaging  under  certain  conditions.  Generating  images 
under  these  conditions  allows  a  comparison  between  the  two  algorithms.  The  results 
shown  are  for  images  that  are  512  x  512  pixels;  used  in  order  to  increase  the  resolution 
of  the  results.  In  Figure  6.6,  at  a  separation  of  2  pixels,  the  polarimeter  deconvolution 
algorithm  has  resolved  the  two  bars  using  the  Sparrow  [16]  criteria  of  a  flat  top  between 
the  two  bars.  The  RL  algorithm  is  able  to  resolve  the  two  bars  at  a  minimum  pixel 
separation  of  4  pixels.  Thus  under  these  lighting  conditions  the  polarimeter  algorithm 
achieves  twice  the  resolution  when  compared  to  the  RL  algorithm. 

Reducing  the  photons  used  for  the  cross  sections  in  Figure  6.6  by  10%  results 
in  a  corresponding  reduction  in  SNR  by  10%  and  results  in  the  set  of  cross  sections 
shown  in  Figure  6.7.  When  the  light  in  the  scenario  is  reduced,  both  algorithms 
require  a  larger  pixel  separation  in  order  to  resolve  the  bar  target.  In  this  case,  the 
RL  algorithm  is  unable  to  resolve  the  bars  until  the  bars,  that  are  only  three  pixels 
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Figure  6.6:  Cross-sectional  views  of  deconvolved  objects  at  an  SNR  of  0.71. 


in  width,  are  separated  by  8  pixels.  The  polarimeter  algorithm  is  able  to  resolve  the 
bars  at  a  bar  separation  of  only  4  pixels. 

As  the  CRLB  results  indicate,  when  lighting  is  sufficiently  high  and/or  as  the 
pixel  separation  is  increased,  the  ability  to  resolve  the  two  bars  becomes  similar. 
Adjusting  the  light  to  achieve  a  SNR  of  6  results  in  the  cross  sectional  views  shown  in 
Figure  6.8.  Both  algorithms  can  resolve  the  bars  at  a  minimum  separation  of  1  pixel. 
The  polarimeter  algorithm  still  has  a  better  fidelity  score  at  the  smaller  separations 
due  to  the  deeper  valley  between  the  peaks.  The  scores  become  nearly  the  same  when 
the  bars  reach  a  separation  of  5  pixels.  However,  the  polarimeter  fidelity  score  is 
better  than  the  RL  algorithm  due  to  the  flatter  trough  between  the  two  bars  and  the 
greater  hairing  that  occurs  at  the  bottom  of  the  RL  bars. 
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Figure  6.7:  Cross-sectional  views  of  deconvolved  objects  at  an  SNR  of  0.65. 


6.1.2  Known  psf  with  laboratory  data.  Using  laboratory  images  produced 
from  the  setup  described  in  Section  4.3,  both  algorithms  are  allowed  to  iterate  un¬ 
til  the  termination  criteria  is  reached.  The  termination  criteria  is  determined  by  the 
amount  of  noise  in  the  images.  Once  the  mean  square  error  between  iterations  reaches 
the  noise  level  present  in  the  images  the  algorithm  terminates.  Both  algorithms  are 
able  to  deconvolve  the  object  of  interest  when  provided  the  estimate  of  the  defocus 
aberration.  The  estimates  of  the  actual  defocus  aberration  used  to  produce  the  results 
are  generated  by  the  FAD  algorithm  of  Section  4.2.  Based  on  the  results  and  the  pre¬ 
vious  discussion  of  the  bounds  in  the  previous  section  this  indicates  that  the  lighting 
is  sufficient  for  both  algorithms  to  have  similar  performance  at  the  lower  aberrations. 
As  the  defocus  aberration  is  increased,  the  polarimeter  performance  increases  corre- 
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Figure  6.8:  Cross-sectional  views  of  deconvolved  objects  at  an  SNR  of  6. 

spondingly  over  the  RL  algorithm.  For  the  laboratory  images,  the  difference  between 
the  two  algorithms  is  seen  in  the  differences  between  the  peaks  and  valleys  of  the 
object  estimates. 


RL  deconvolved  object 


50  100  150  200  250 


Pol  deconvolved  object 
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RL  deconvolved  object 


Pol  deconvolved  object 


Figure  6.9:  Deconvolution  -  Position  1. 
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In  Figure  6.9  the  RL  algorithm’s  first  valley  is  39%  of  the  first  peak  and  the 
resulting  area  between  the  bars  in  the  deconvolved  object  is  brighter.  The  bars  are 
less  distinguishable.  The  polarimeter  algorithm’s  first  valley  is  24%  of  the  first  peak 
and  results  in  a  much  darker  area  between  the  bars  making  them  easier  to  distinguish 
their  extent. 
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Figure  6.10:  Deconvolution  -  Position  5. 


In  Figure  6.10  the  RL  algorithm’s  first  valley  is  38%  of  the  first  peak  and  the 
polarimeter  algorithm’s  first  valley  is  30%  of  the  first  peak. 

The  target  is  nearly  in  focus  in  Figure  6.11  at  Position  7  and  both  algorithms 
have  very  similar  performance.  This  is  consistent  with  the  lower  aberration  and  ample 
lighting.  The  laboratory  results  confirm  the  fidelity  predictions  of  the  simulation 
results  of  Section  6.1.  The  lower  defocus  aberration  results  of  both  algorithms  would 
be  scored  nearly  the  same  but  the  polarimeter  algorithm  would  receive  a  better  score 
as  the  defocus  aberration  is  increased.  The  difference  between  the  two  algorithms  is 
more  significant  when  considering  the  results  of  their  use  in  the  multiframe  algorithms 
in  the  following  sections. 
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Figure  6.11:  Deconvolution  -  Position  7. 


6.2  Polarimeter  Multiframe  Blind  Deconvolution 

This  section  presents  the  results  of  the  extension  of  both  the  RL  algorithm  to 
the  MFBD  algorithm  and  the  polarimeter  algorithm  to  the  PMFBD  algorithm.  The 
simulation  uses  three  separate  primary  channel  images,  Figure  6.12,  for  testing  the 
MFBD  algorithm.  The  six  corresponding  polarimeter  images,  3  primary  channels  and 
3  polarization  sensitive,  are  shown  in  Figure  6.13. 


(a)  Image  1.  (b)  Image  2.  (c)  Image  3. 

Figure  6.12:  Images  used  for  MFBD  algorithm  testing. 


The  three  channels  vary  by  different  amounts  of  defocus  aberration.  Since  the 
images  are  simulated  there  is  no  need  to  register  the  independent  images  resulting 
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(d)  Pol  Image  1.  (e)  Pol  Image  2.  (f)  Pol  Image  3. 

Figure  6.13:  Images  used  for  PMFBD  algorithm  testing. 


in  zero  registration  error.  Image  1  has  an  aberration  of  150  times  the  atmospheric 
defocus  covariance  with  a  D/R0  =  1.  Image  2  has  an  aberration  of  90  times  the 
atmospheric  defocus  covariance  with  a  D/R0  =  1.  Image  3  has  an  aberration  of  40 
times  the  atmospheric  defocus  covariance  with  a  D/R0  =  1.  The  MFBD  images  have 
twice  the  SNR  as  the  primary  images  for  the  PMFBD  simulation  to  account  for  the 
theoretical  beamsplitter.  This  is  visible  in  the  different  amounts  of  noise  visible  in 
the  images. 

The  initial  guess  for  both  algorithms  is  a  diffraction  limited  OTF  for  all  three 
objects.  Both  algorithms  then  follow  a  similar  breakdown  in  the  internal  iteration. 
First  they  estimate  the  new  objects  from  the  old  objects  using  the  old  estimates  of 
the  psfs.  Then  they  update  the  three  psfs  estimates  using  the  old  object  estimates. 

In  order  to  compare  performance  all  internal  loops  arc  limited  to  a  fixed  number 
of  iterations.  For  instance,  each  object  estimate  update  is  limited  to  20  iterations. 
The  psf  updates  are  limited  to  5  iterations  each.  The  outside  loop  is  allowed  to 
converge  to  a  fixed  mean  square  error  between  executions  of  the  outside  loop.  This 
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(a)  Estimated  object. 


Pixel  location 

(b)  Object  cross-section. 


Figure  6.14:  MFBD  estimated  object. 


- Actual 

- Estimate 


Figure  6.15:  MFBD  estimated  OTFs. 


predetermined  value  is  based  on  the  noise  level  of  the  images.  For  the  simulations, 
the  standard  deviation  of  the  noise  is  approximately  0.6.  Using  this  value,  the  mean 
square  error  (MSE)  between  outside  loops  is  allowed  to  converge  to  a  value  of  0.7. 
The  overall  effect  is  that  the  MFBD  and  PMFBD  algorithms  progress  in  lock  step 
based  on  number  of  object  and  psf  updates.  The  MFBD  algorithm  usually  requires 
a  significantly  higher  number  of  outer-loop  iterations  to  converge  to  the  same  MSE 
error  between  outer  loops. 

The  MFBD  results  are  shown  in  Figures  6.14  and  6.15.  Even  with  a  target  bar 
separation  of  8  pixels  the  MFBD  algorithm  has  a  hard  time  deconvolving  the  two 
bars.  There  is  a  visible  valley  between  the  two  sides  of  the  object.  The  three  psf 
estimates  are  Fourier  transformed  into  their  OTFs  in  order  to  visually  compare  them 
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Pixel  location 

(a)  Estimated  object.  (b)  Object  cross-section. 

Figure  6.16:  PMFBD  estimated  object. 


- Actual 

- Estimate 


Figure  6.17:  PMFBD  estimated  OTFs. 

to  the  actual  OTFs  used  to  create  the  images.  The  results  seen  in  Figure  6.15  show 
that  the  MFBD  algorithm  is  able  to  estimate  OTFs  that  are  closer  to  the  actual  OTFs 
than  the  original  guess  at  the  start  of  the  algorithm.  However,  the  error  between  the 
estimates  and  the  actual  OTFs  is  still  significant  when  the  algorithm  reaches  the  MSE 
termination  criteria. 

The  results  of  the  PMFBD  are  shown  in  Figures  6.16  and  6.17.  The  simulated 
bar  target  is  clearly  visible  in  Figure  6.16(a).  The  cross  section  visible  in  Figure 
6.16(b)  shows  a  clear  separation  in  the  bars.  These  results  are  a  clear  improvement 
over  the  MFBD  algorithm  with  the  same  input  images  even  with  the  lower  light  levels 
present  the  polarimeter  images.  The  significant  improvement  is  due  to  the  much  better 
estimation  of  the  psfs.  Looking  at  the  estimated  OTFs  in  Figure  6.17,  there  is  a  very 
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close  match  between  the  estimated  OTF  and  the  actual  OTF.  The  polarization  data 
allows  a  much  better  estimate  of  the  psf  to  be  made. 

6.2.1  Convergence  Rates.  With  both  algorithms  allowed  a  fixed  number 
of  internal  iterations  per  outer  loop,  the  convergence  to  the  same  MSE  termination 
condition  is  presented.  The  results  for  the  previous  section’s  same  three  images  are 
presented  in  Figure  6.18.  The  algorithms  are  run  consecutively  for  bar  separations 
from  1  to  8  pixels.  The  same  bar  width  of  3  pixels  is  used.  The  PMFBD  algorithm 
converges  to  a  much  better  solution  based  on  the  fidelity  metric  within  half  the  outer 
loop  iterations  of  the  MFBD  algorithm.  The  MFBD  algorithm  fails  to  make  significant 
improvements  in  object  estimation  by  the  time  the  termination  criteria  is  reached. 


(a)  PMFBD  convergence.  (b)  MFBD  convergence. 

Figure  6.18:  Algorithm  fidelity  versus  iterations. 

6.2.2  Laboratory  Results.  The  three  laboratory  images  at  Position  3,  Po¬ 
sition  4,  and  Position  5  are  given  to  both  the  MFBD  and  the  PMFBD  algorithms. 
The  resulting  estimate  of  the  object  is  shown  in  Figure  6.19.  The  significant  differ¬ 
ence  in  object  estimation  between  the  two  algorithms  is  in  the  peaks  and  valleys  of 
the  object  estimated.  The  MFBD  object  valleys  are  approximately  30%  of  the  peak 
values.  The  PMFBD  valleys  almost  go  to  zero  given  the  same  images.  The  valleys 
are  approximately  1  /  13th  of  the  peak  values.  The  actual  peak  values  are  very  similar. 
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(a)  MFBD  Estimate.  (b)  PMFBD  Estimate. 


(c)  MFBD  cross-section. 


(d)  PMFBD  cross-section. 


Figure  6.19:  Laboratory  multiframe  estimated  objects  and  cross-sections. 


The  most  significant  difference  between  the  two  algorithms  is  in  the  ability  to 
estimate  the  actual  OTFs.  The  results  of  the  OTF  estimates  are  shown  in  Figure 
6.20.  The  top  set  of  estimates  is  produced  by  the  MFBD  algorithm  and  the  bottom 
set  is  produced  by  the  PMFBD  algorithm.  The  actual  OTF  used  for  comparison  is  the 
OTF  estimated  by  the  FAD  algorithm  for  the  laboratory  images.  With  similar  results 
to  simulation,  the  MFBD  algorithm  fails  to  significantly  improve  the  estimates  of  the 
OTFs  for  the  three  objects.  The  PMFBD  algorithm  clearly  improves  the  estimation 
of  the  OTFs  for  the  three  objects. 

The  results  in  this  chapter  show  that  polarization  data,  when  available,  improves 
the  ability  to  deconvolve  objects  using  the  algorithms  developed  for  polarimeter  de- 
convolution  and  polarimeter  blind  deconvolution. 
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- Actual 

- Estimate 


(a)  MFBD  estimated  OTF  for  (b)  MFBD  estimated  OTF  for  (c)  MFBD  estimated  OTF  for 
object  1.  object  2.  object  3. 


(d)  PMFBD  estimated  OTF  (e)  PMFBD  estimated  OTF  for  (f)  PMFBD  estimated  OTF  for 
for  object  1.  object  2.  object  3. 


Figure  6.20:  Laboratory  multiframe  estimated  OTFs. 


Ill 


VII.  Conclusions 


The  previous  chapters  detail  the  objectives  of  this  research  and  their  successful 
achievement.  The  first  objective  achieved  is  detailed  in  Chapter  Iff  which  shows 
the  development  of  the  spatial  resolution  CRLB  for  a  two-channel  polarimeter  incor¬ 
porating  atmospheric  effects.  The  second  objective  is  achieved  with  the  production 
of  calibrated  data  detailed  in  Chapter  IV.  The  calibrated  data  is  used  for  testing  of 
the  algorithms  developed  in  Chapter  V  for  the  third  objective. 

The  spatial  resolution  CRLB  for  a  two-channel  polarimeter  show  improved  spa¬ 
tial  resolution  under  certain  lighting  conditions.  With  sufficient  light  both  systems 
have  the  same  performance.  However,  below  a  certain  point  as  the  light  levels  de¬ 
crease,  the  polarimeter’s  spatial  resolution  improves  over  standard  imaging  as  the 
degree  of  polarization  of  the  object  increases.  The  limit  in  spatial  resolution  achieved 
by  either  imaging  system  is  determined  by  the  light  level.  The  conditions  in  which  the 
polarimeter  outperforms  standard  imaging  is  shown  in  the  comparison  of  the  CRLBs 
for  both  systems.  Thus  the  first  objective  is  achieved. 

The  second  objective  is  achieved  by  producing  a  set  of  calibrated  data  for  testing 
purposes.  An  optics  bench  is  set  up  to  provide  an  initial  set  of  real  images  with  con¬ 
trolled  focus  errors.  The  FAD  algorithm,  presented  in  Section  4.2,  is  a  new  approach 
for  detection  of  focus  error  using  a  single  image.  This  is  an  improvement  over  histor¬ 
ical  methods  that  require  complicated  calibrations  or  the  diversion  of  light  from  the 
primary  imaging  channel.  The  ability  to  detect  algorithmically  the  actual  focus  error 
in  an  image  allows  for  more  accurate  test  images  from  the  laboratory  setup.  This 
reduces  the  error  in  location  of  the  actual  focal  point  of  the  laboratory  setup  that 
would  otherwise  result  from  the  use  of  real  (non-thin)  lenses,  a  compound  lens,  and 
locations  of  other  components.  The  laboratory  images  were  used  to  develop  Matlab® 
images  with  similar  statistics  but  with  a  broader  range  of  lighting  conditions.  Both 
the  laboratory  images  and  the  simulated  images  were  useful  in  testing  the  algorithms 
developed. 
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The  third  objective  is  achieved  with  the  successful  development  of  the  algorithms 
in  Chapter  V  that  make  use  of  the  polarization  information.  Initially,  an  algorithm  is 
developed  for  polarimeter  deconvolution  with  a  known  aberration.  The  incorporation 
of  the  polarization  information  yields  an  improvement  in  deconvolution.  One  can 
see  from  the  cross-sectional  views  of  Figure  6.6  and  6.7  that  the  improved  resolution 
based  on  the  bound  calculations  of  Chapter  III  is  realized.  At  a  SNR  of  0.71,  Figure 
6.6  shows  that  the  polarimeter  can  resolve  the  two  bars  at  a  minimum  separation 
of  two  pixels  while  the  standard  imaging  system  resolves  the  bars  with  a  minimum 
of  four  pixels.  Thus  the  polarimeter  achieves  half  the  spatial  resolution  compared 
to  the  standard  imaging  system.  Figure  6.7  shows  that  a  further  reduction  in  the 
number  of  photons  results  in  the  failure  of  the  RL  algorithm  to  resolve  the  bars  but 
the  polarimeter  algorithm  resolves  the  bars  at  a  separation  of  four  pixels.  As  the 
pixel  separation  and/or  lighting  increases,  the  ability  to  resolve  the  two  bars  becomes 
nearly  identical.  These  results  conform  to  the  CRLB  comparisons  between  the  two 
systems. 

An  additional  benefit  of  the  polarimeter  deconvolution  is  seen  in  the  fidelity 
metric  results  shown  in  Figures  6.2,  6.3,  and  6.4.  The  polarimeter  deconvolution  con¬ 
sistently  has  a  higher  deconvolved  object  fidelity  than  the  RL  algorithm  when  they 
are  operating  on  the  same  images.  The  reason  for  this  is  the  polarimeter  deconvolu¬ 
tion’s  ability  to  estimate  bars  with  steeper  sides  and  lower  valleys  when  compared  to 
the  RL  algorithm  results. 

The  polarization  deconvolution  algorithm  is  then  extended  under  the  assump¬ 
tion  of  an  unknown  aberration.  Leveraging  the  results  of  Schulz’  MFBD  algorithm 
results  in  the  PMFBD  algorithm.  The  PMFBD  algorithm  is  a  GEM  algorithm  that 
alternates  between  estimation  of  the  object  and  estimation  of  the  psf.  The  PMFBD 
algorithm  achieves  significant  improvement  in  the  rate  of  convergence  over  the  RL 
based  MFBD  algorithm,  Section  6.2.1.  The  most  significant  improvement  over  the 
MFBD  algorithm  is  in  the  estimation  of  the  OTFs.  The  OTFs  estimated  by  the 
PMFBD  algorithm,  Figure  6.17,  are  very  close  to  the  OTFs  used  to  create  the  test 
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images.  However,  the  OTFs  estimated  by  the  MFBD  algorithm  are  significantly  dif¬ 
ferent,  as  in  Figure  6.15.  The  addition  of  the  polarization  information  in  the  PMFBD 
algorithm  produces  superior  results. 

All  of  these  results  put  together  show  that  the  use  of  polarization  data  in  electro- 
optical  imaging  systems  improves  the  ability  to  deconvolve  objects  of  interest.  This 
research  goal  is  clearly  achieved. 

Recommended  Research 

Recommended  follow-ons  to  this  research  are: 

•  To  extend  the  polarimeter  deconvolution  algorithms  with  additional  channels. 
With  only  two  channels  currently  used,  one  cannot  distinguish  between  un¬ 
polarized  light  and  light  polarized  at  45  degrees  with  respect  to  the  polarized 
channel.  Additional  polarimeter  channels  would  allow  estimation  of  the  degree 
of  polarization  present  in  the  scene. 

•  To  implement  the  FAD  algorithm  using  a  Weiner  filter  to  deconvolve  the  object 
instead  of  the  RL  algorithm,  as  discussed  in  Section  4.4.3. 
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